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Abstract 


A new high-resolution and genuinely multidimensional numerical method for solving 
conservation laws is being developed. It was designed to avoid the limitations of the 
traditional methods, and was built from ground zero with extensive physics considerations. 
Nevertheless, its foundation is mathematically simple enough that one can build from it a 
coherent, robust, efficient and accurate numerical framework. 

Two basic beliefs that set the new method apart from the established methods are 
at the core of its development. The first belief is that, in order to capture physics more 
efficiently and realistically, the modeling focus should be placed on the original integral 
form of the physical conservation laws, rather than the differential form. The latter form 
follows from the integral form under the additional assumption that the physical solution 
is smooth, an assumption that is difficult to realize numerically in a region of rapid change, 
such as a boundary layer or a shock. The second belief is that, with proper modeling of the 
integral and differential forms themselves, the resulting numerical solution should auto- 
matically be consistent with the properties derived from the integral and differential forms, 
e.g., the jump conditions across a shock and the properties of characteristics. Therefore a 
much simpler and more robust method can be developed by not using the above derived 
properties explicitly. 

Specifically, to capture physics as fully as possible, the method requires that: (i) space 
and time be unified and treated as a single entity; (ii) both local and global flux conser- 
vation in space and time be enforced; and (iii) a multidimensional scheme be constructed 
without using the dimensional-splitting approach, such that multidimensional effects and 
source terms (which are scalars) can be modeled more realistically. 

To simplify mathematics and broaden its applicability as much as possible, the method 
attempts to use the simplest logical structures and approximation techniques. Specifically, 
(i) it uses a staggered space-time mesh such that flux at any interface separating two con- 
servation elements can be evaluated internally in a simpler and more consistent manner, 
without using a separate flux model; (ii) it does not use many well-established techniques 
such as Riemann solvers, flux splittings and monotonicity constraints such that the limi- 
tations and complications associated with them can be avoided; and (iii) it does not use 
special techniques that are not applicable to more general problems. 

Furthermore, triangles in 2D space and tetrahedrons in 3D space are used as the basic 
building blocks of the spatial meshes, such that the method (i) can be used to construct 
2D and 3D non-dissipative schemes in a natural manner; and (ii) is compatible with the 
simplest unstructured meshes. 

Note that while numerical dissipation is required for shock capturing, it may also result 
in annihilation of small disturbances such as sound waves and, in the case of flow with a 
large Reynolds number, may overwhelm physical dissipation. To overcome this difficulty, 
two different and mutually complementary types of adjustable numerical dissipation are 
introduced in the present development. 


N AS A/TM— 1 998-208843 


1 



1. Introduction 


Since its inception in 1991 [1], the space-time conservation element and solution ele- 
ment method [1-32] has been used to obtain highly accurate numerical solutions for flow 
problems involving shocks, rarefaction waves, acoustic waves, vortices, ZND detonation 
waves, shock/acoustic waves/ vortices interactions, dam- break and hydraulic jump. This 
article is the first of a series of papers that will provide a systematic and up-to-date descrip- 
tion of this new method (hereafter it may be referred to abbreviatedly as the space-time 
CE/SE method or simply as the CE/SE method). To answer frequently-asked questions 
and clarify possible misconceptions, we shall begin this paper with (i) an overall view of the 
CE/SE method and its capabilities, and (ii) an extensive comparison of the basic concepts 
used by the CE/SE method with those used by other methods. 

Currently, the field of computational fluid dynamics (CFD) represents a diverse col- 
lection of numerical methods, with each of them having its own limitations. Generally 
speaking, these methods were originally introduced to solve special classes of flow prob- 
lems. Development of the CE/SE method is motivated by a desire to build a brand new, 
more general and coherent numerical framework that avoids the limitations of the tradi- 
tional methods. 

The new method is built on a set of design principles given in [2]. They include: (i) 
To enforce both local and global flux conservation in space and time, with flux evaluation 
at an interface being an integral part of the solution procedure and requiring no interpo- 
lation or extrapolation; (ii) To unify space and time and treat them as a single entity; 
(iii) To consider mesh values of dependent variables and their derivatives as independent 
variables, to be solved for simultaneously; (iv) To use only local discrete variables rather 
than global variables like the expansion coefficients used in spectral methods; (v) To de- 
fine conservation elements and solution elements such that the simplest stencil will result; 
(vi) To require that, as much as possible, a numerical analogue be constructed so as to 
share with the corresponding physical equations the same space-time invariant properties, 
such that numerical dissipation can be minimized [5,10,24]; (vii) To exclude the use of 
characieristics-based techniques (such as Riemann solvers ); and (viii) To avoid the use of 
ad hoc techniques as much as possible. 

Moreover, the development of the CE/SE method is also guided by two basic beliefs 
that set it apart from the established methods. The first belief is that, in order to cap- 
ture physics more efficiently and realistically, the modeling focus should be placed on the 
original integral form of the physical conservation laws, rather than the differential form. 
The latter form follows from the integral form under the additional assumption that the 
physical solution is smooth, an assumption that is difficult to realize numerically in a re- 
gion of rapid change , such as a boundary layer or a shock . The second belief is that, with 
proper modeling of the integral and differential forms themselves, the resulting numerical 
solution should automatically be consistent with the properties derived from the integral 
and differential forms, e.g., the jump conditions across a shock and the properties of char- 
acteristics. In other words, a much simpler and more robust method can be developed by 
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not using the above derived properties explicitly. 

With the exception of the Navier-Stokes solver, all the ID schemes described in [2] 
have been extended to become their 2D counterparts [9-11,14]. A more complete account 
of these new 2D schemes and their applications will be given in this and the following 
papers [3,4]. It will be shown in Sec. 3 that the spatial meshes used in these schemes 
are built from triangles — in such a manner that the resulting meshes are completely dif- 
ferent from those used in the finite element method. As a result, these schemes are (i) 
compatible with the simplest unstructured meshes [31], and (ii) constructed without using 
the dimensional-splitting approach, i.e., without applying a ID scheme in each coordinate 
direction. The dimensional-splitting approach is widely used in the construction of multidi- 
mensional upwind schemes. Unfortunately, this approach is flawed in several respects [33]. 
In particular, because a source term is not aligned with a special direction, it is difficult to 
imagine how this dimensional-splitting approach, in a logically consistent manner, can be 
used to solve a multidimensional problem involving source terms, such as those modeling 
chemical energy release. 

Moreover, as will be shown shortly, because the CE/SE 2D schemes share with their 
ID versions the same design principles, not only is the extension to 2D a straightforward 
matter, each of the new 2D schemes also shares with its ID version virtually identical 
fundamental characteristics. 

At this juncture, note that monotonicity conditions are not observed by general flow 
fields, e.g., those involving ZND detonation waves [21]. As a result, techniques involving 
monotonicity constraints are not used in the present development. 

To give the reader, in advance, a concrete example that demonstrate the validity 
of the two basic beliefs referred to earlier, a self-contained Fortran program is fisted in 
Appendix A. It is a CE/SE solver [23] for an extended Sod’s shock tube problem that 
is the original Sod’s problem [38] with the additional complication of imposing a non- 
reflecting boundary condition at each end of the computational domain. Note that the 
flow under consideration contains discontinuities and, relative to the computational frame, 
is subsonic throughout. It is well known that implementing a non-reflecting boundary 
condition for a subsonic flow is much more difficult than doing the same for a supersonic 
flow. This difficulty is further exacerbated by the fact that the traditional non-reflecting 
boundary conditions, e.g., the characteristic, the radiation (asymptotic), the buffer-zone, 
and the absorbing boundary conditions [39-44] are all based on an assumption that is not 
valid for the present case, i.e., that the flow is continuous. In spite of the fact that solving 
the present extended Sod’s problem is substantially more difficult than the original Sod s 
problem, the main loop in the program fisted herein contains only 39 Fortran statements. 
Not only is it very small in size, this program has a very simple logical structure. With 
the exception of a single “if” statement used to identify the time levels at which the 
non-reflecting boundary conditions must be imposed, it contains no conditional Fortran 
statements or functions such as “if”, “amax”, or “amin” that are often used in programs 
implementing high-resolution upwind methods. The small size of the fisted program reflects 
the simplicity of the techniques employed by the CE/SE method to capture shock waves. 
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It also results from the fact that the non-reflecting boundary conditions used in the present 
solver are the simple extrapolation conditions Eqs. (2.66) and (2.67) given in Sec. 2. They 
are much simpler than the traditional non-reflecting boundary conditions. On the other 
hand, the absence of Fortran conditional statements is a result of avoiding the use of 
ad hoc techniques. In spite of its small size and simple logical structure, according to 
the numerical results generated by the listed program (presented here as Figs. l(a)-(c), 
with the numerical results and the exact solutions denoted by triangles and solid lines, 
respectively; see also [23]), the present solver is capable of generating nearly perfect non- 
reflecting solutions using the same time-step size from t = 0. Note that, at t = 10, all the 
waves have exited the computational domain, i.e., the exact solution is constant within 
it. The theoretical values of density, velocity, and pressure are approximately 0.4262000, 
0.9277462 and 0.3030000, respectively. The maximum magnitudes of the errors in the 
numerically computed values of density, velocity, and pressure at t = 10 are approximately 
0.0004, 0.0007, and 0.0004, respectively. 

Note that Eqs. (2.66) and (2.67) represent only one of many sets of simple and ro- 
bust non-reflecting boundary conditions developed especially for the CE/SE method [23]. 
Behind this development is a radical new concept based entirely on an assumption about 
the space-time flux distribution in the neighborhood of a spatial boundary. As a result, 
implementation of these CE/SE non-reflecting boundary conditions does not require the 
use of characteristics-based techniques. 

To further demonstrate the nontraditional nature of the CE/SE method, the numeri- 
cal results generated using the steady-state non-reflecting boundary conditions that were 
introduced and rigorously justified in [23] will also be presented here. Consider an alter- 
nate CE/SE solver that differs from the above CE/SE solver only in the fact that the 
steady-state boundary conditions Eq. (2.68) given in Sec. 2 are now taking the place of 
Eqs. (2.66) and (2.67). At t = 0.2, the waves generated in the interior of the computational 
domain have not yet reached the boundaries. In this case, with the given initial conditions 
(i.e., two different uniform states separated by a discontinuity located at the dead center of 
the domain), each of the above two solvers yield the same uniform solution in the vincinity 
of the right or left boundary. As a result, at t = 0.2, the numerical results generated by 
the alternate solver are identical to those shown in Fig. 1(a). The numerical results of the 
alternate solver at t = 0.4 are shown in Fig. 2(a). It is seen that, by this time, the shock 
wave has passed cleanly through the right boundary. There is good agreement between 
the numerical solution and the exact solution everywhere in the interior except for a slight 
disagreement in the vicinity of the right boundary. Note that the right boundary values, 
which do not vary with time, are discontinuous with respect to the neighboring interior 
values. The numerical results at t = 0.6 are shown in Fig. 2(b). As seen from the density 
profile, by this time, the contact discontinuity has also passed through the right boundary. 
Agreement between the numerical solution and the exact solution continue to be good in 
the interior. However, both left and right boundary values are now discontinuous with 
respect to the neighboring interior values. 

Note that several recent applications [13,16,17,26,28] of the CE/SE method to 2D 


NASA/TM — 1998-208843 


4 



aeroacoustics problems reveal that: (i) the trivial nature of implementing CE/SE non- 
reflecting boundary conditions is manifested even for 2D problems; (ii) accuracy of the 
numerical results for nonlinear Euler problems is comparable to that of a 4-6th order 
compact difference scheme, even though nominally the CE/SE solver used is only of 2nd- 
order accuracy; and (iii) most importantly, the CE/SE method is capable of accurately 
modeling both small disturbances and strong shocks, and thus provides a unique tool 
for solving flow problems where the interactions between sound waves and shocks are 
important, such as the noise field around a supersonic over- and under-expanded jet. The 
fact listed in item (i) is more fundamental in nature, and will be further discussed in a 
separate paper. The following comments pertain to items (ii) and (iii): 

(a) Assuming the same order of accuracy, generally speaking, the accuracy of a scheme 
that enforces the space-time flux- conservation property is higher than that of a scheme 
that does not. A compact scheme generally does not enforce the flux-conservation 
property of the nonlinear Euler equations. On the contrary, not only is the present 
scheme flux-conserving, its accuracy in nonlinear calculations is enhanced by its sur- 
prisingly small dispersive errors [2,8,13,16,17]. Moreover, the nominal order of accu- 
racy of an Euler solver is determined assuming a linearized form of the Euler equations. 
Thus its significance with respect to a highly nonlinear solution of the Euler equations 
may be questionable. 

(b) while numerical dissipation is required for shock resolution, it may also result in 
annih ilation of small disturbances such as sound waves. Thus, a solver that can handle 
both small disturbances and strong shocks must be able to overcome this difficulty. 
It will be explained shortly that the CE/SE method is intrinsically endowed with this 
capability. On the other hand, a high-resolution upwind scheme that focuses only on 
shock resolution may introduce too much numerical dissipation [45]. 

Next we shall review briefly the inviscid version of the a-p scheme described in [2]. In 
addition to providing a historical perspective, the review will remove, once and for all, any 
lingering doubt from the reader’s mind that the CE/SE method indeed differs substantially 
in both concept and methodology from the well-established methods. In particular, it will 
give in advance answers to questions such as: (i) is there any difference between the space- 
time elements used here and those used in the finite element method? and (ii) what are 
the key differences between the CE/SE method and other finite volume methods? 

To proceed, consider an initial- value problem involving the PDE 


du du 

di +a dx 


= 0 


( 1 . 1 ) 


where the convection speed a is a constant. The exact solution to any such problem 
has three fundamental properties: (i) it does not dissipate with time; (ii) its value at a 
spatial point at a later time has a finite domain of dependence (a point) at an earlier 
time; and (iii) it is completely determined by the initial data at a given time. Ideally, a 
numerical solution for Eq. (1.1) should also possess the same three properties. Because 
(i) a solution of a dissipative numerical scheme will dissipate with time, (ii) the value of a 
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solution of an implicit scheme at any point ( x,t ) is dependent on all initial data, and all 
the boundary data up to the time t , and (iii) the unique determination of a solution by a 
scheme involving more than two time levels requires the specification of the data at at least 
the first two time levels, an ideal solver must be a two-level, explicit, and non- dissipative 
(i.e., neutrally stable) scheme. In 1991, the first solver known to the authors that satisfies 
the above conditions was reported in [1]. Because this new solver models Eq. (1.1) which 
is characterized by the parameter a, it is referred to as the a scheme. The a scheme is 
non-dissipative if the Courant number is less than unity. 

At this juncture, the reader may wonder what the merit is of constructing a neutrally 
stable scheme. After all, it is well known that its nonlinear extensions generally are unsta- 
ble. To address this question, the significance of constructing such a scheme and the critical 
role it plays in the development of the CE/SE method will be discussed immediately. 

To proceed, note that there are several explicit and implicit extensions [2,12,25] of the 
a scheme which are solvers for 


du du d 2 u 

m +a ax~ t ‘a^ = 0 


( 1 . 2 ) 


Here the viscosity coefficient p(> 0) is a constant. Because Eq. (1.2) is characterized by 
the parameters a and p, these extensions are referred to as either the explicit a-p schemes 
or the implicit a-p schemes. Each of these schemes reduces to the non-dissipative a scheme 
when p = 0. As a result, each of them has the property that the numerical dissipation of 
its solutions approaches zero as the physical dissipation approaches zero. 

The above property is important because of the following observation: with a few 
exceptions, the numerical solution of a time-marching problem generally is contaminated 
by numerical dissipation. For a nearly inviscid problem, e.g., flow at a large Reynolds 
number, numerical dissipation may overwhelm physical dissipation and cause a complete 
distortion of the solution. To avoid such a difficulty, ideally a CE/SE solver for Eq. (1.2) 
or for the Navier-Stokes equations should possess the above special property. Obviously 
the development of such a solver must be preceded by that of a neutrally stable solver of 
Eq. (1.1). 

The problem of physical dissipation being overwhelmed by numerical dissipation does 
not exist for a pure convection problem. However, as explained in the earlier discussion 
about the delicate nature of simulating small disturbances in the presence of shocks, nu- 
merical dissipation must still be handled carefully in this case. 

Note that numerical dissipation traditionally is adjusted by varying the magnitude of 
added artificial dissipation terms. However, after being stripped of these added artificial 
dissipation terms, almost every traditional scheme (such as the Lax-Wendroff scheme ) is 
still not free from inherent numerical dissipation. Hence, numerical dissipation generally 
cannot be avoided completely using the traditional approach. 

This completes the discussion about the roles of non-dissipative schemes in the current 
development. To proceed further, the construction of the ID a scheme will be described 
briefly. 
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Let Xi = x , and x^ = t be considered as the coordinates of a two-dimensional Eu- 
clidean space E 2 ■ By using Gauss’ divergence theorem in the space- time E 2 , it can be 
shown that Eq. (1.1) is the differential form of the integral conservation law 

([ h • ds = 0 (1*3) 

Js(V) 

As depicted in Fig. 3, here (i) S(V) is the boundary of an arbitrary space-time region V 
in E 2 ; (ii) h = ( au , u) is a current density vector in £ 2 ; and (iii) ds = dan with da and 
n, respectively, being the area and the outward unit normal of a surface element on S(V). 
Note that (i) h • ds is the space-time flux of h leaving the region V through the surface 
element ds, and (ii) all mathematical operations can be carried out as though E 2 were an 
ordinary two-dimensional Euclidean space. 

Let fl denote the set of all mesh points (j, n) in E 2 (dots in Fig. 4(a)) with n being 
a half or whole integer, and ( j — n) being a half integer. For each (j, n) E Q, let the 
solution element SE(j,n) be the interior of the space-time region bounded by a dashed 
curve depicted in Fig. 4(b). It includes a horizontal line segment, a vertical line segment, 
and their immediate neighborhood. For the discussions given in this paper, the exact size of 
this neighborhood does not matter. However, in case the conservation law Eq. (1.3) takes 
a more complicated form in which the right side is a volume integral involving a source 
term, the SEs must fill the entire computational domain such that the volume integral can 
be modeled properly [21,22]. A SE that fulfills this requirement is depicted in Fig. 4(c). 

For any (x,t) E SE (j,n), let u{x,t) and h(x,t), respectively, be approximated by 
u*{x,t\j,n) and h*(x,t\j,n) which we shall define shortly. Let 

u*(x,t ',j,n) = u” + {u x Yj(x — Xj) + { u t)j (t ~~ t ) (1-4) 

where (i) uj, («*)", and (u*)” are constants in SE(j,n), and (ii) ( xj,t n ) are the coordinates 
of the mesh point (j, n). 

We shall require that u = u*{x,t;j,n ) satisfy Eq. (1.1) within SE(j,n). As a result, 

(u t )? = -a(tt*)7 (1-5) 


Combining Eqs. (1.4) and (1.5), one has 

u*{x,t-,j,n) = u] + (u x )][{x-x j )-a(t-t n )}, (x,t) E SE(j,n) (1.6) 

As a result, there are two independent marching variables u ” and ( tii )J associated with 
each (j, n) E fi. Furthermore, because h = ( au,u ), we define 

h*(x,t ; j,n) = ( au*(x,t;j,n ), u*(x,t ; j,n)) (1-7) 
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Let E 2 be divided into non-overlapping rectangular regions (see Fig. 4(a)) referred to 
as conservation elements (CEs). As depicted in Figs. 4(d) and 4(e), the CE with its top- 
right (top-left) vertex being the mesh point ( j,n ) G Q is denoted by CE_(j,n) (CE + (j,n.)). 
The discrete approximation of Eq. (1.3) is then 

( f h* -ds = 0 (1.8) 

JS(CE±(j,n)) 


for all (j,n) 6 fl. At each (j, n) G fi, Eq. (1.8) provides the two conditions needed to 
solve its two independent marching variables. In the following, the manner in which the 
integrals in Eq. (1.8) should be evaluated will be explained by considering the case that 
involves CE_(j,n). 

According to Fig. 4(d), S(CE-(j, n)), i.e., the boundary of CE_(jf,n), is formed by 
four line segments. Among them, AB and AZ? lie within SE(ji,n). As a result, the 
flux leaving CE_(j,n) through these two line segments will be evaluated using Eqs. (1.6) 
and (1.7) with the assumption that any point (®,<) on them belongs to SE(j,n). On 
the other hand, because CB and CD He within SE(y — 1/2 , n — 1/2), the flux leaving 
CE -(j,n) through them will be evaluated assuming any point ( x,t ) on them belongs to 
SE(j — 1/2, n — 1/2). 

According to Eq. (1.8), the total flux of h* leaving the boundary of any conserva- 
tion element is zero. Because the surface integration over any interface separating two 
neighboring CEs is evaluated using the information from a single SE, obviously the lo- 
cal conservation relation Eq. (1.8) leads to a global flux conservation relation, i.e., the 
total flux of h* leaving the boundary of any space-time region that is the union of any 
combination of CEs will also vanish. 

From the above discussions, it becomes obvious that the space-time element used in 
the Unite element method differs from the current space-time SE and CE in both concept 
and the roles they serve. In particular, the former is not introduced to enforce flux conser- 
vation. In contrast to this, in the CE/SE method, flux conservation transmits information 
between neighboring SEs, and no global smoothness requirements are made on the solu- 
tion to link neighboring SEs. This strategy enables the accurate capturing of traveling 
multidimensional solution discontinuities, e.g., moving multidimensional shock waves. 

Furthermore, the CE/SE method is also fundamentally different from the traditional 
finite-volume methods such as the high-resolution upwind methods [46,47] and the dis- 
continuous Galerkin method [48] in one important respect, i.e., because of the space-time 
staggering nature of its solution elements, the present method has a much simpler and 
consistent procedure to evaluate the dux at an interface. The key features of CE/SE flux- 
evaluation that distinguish it from those of the traditional methods are discussed in the 
following remarks: 

(a) Because an interface separating two neighboring CEs lies within a SE, the flux at this 

interface is evaluated without interpolation or extrapolation. Furthermore, the SE to 
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which a particular interface belongs is determined by a rule that is independent of the 
local numerical solution. In other words, the concept of special upwind treatments 
and the complications that arise from these treatments are entirely f oreig n to the 
CE/SE method. To be more specific, consider the flux at the interface AD depicted 
in Fig. 4(d). It is completely determined by uj and (u x )J, two numerical variables 
associated with the predetermined mesh point (j,n), i.e., point A. 

(b) Flux evaluation is straightforward and it requires only simple integration involving the 
first-order Taylor’s expansion. No comphcated techniques such as the characteristics- 
based techniques are ever needed. 

Finally, we also want to emphasize that the concepts used in the construction of the 
a scheme are fundamentally different from several schemes introduced by Nessyahu and 
Tadmor[49], and Sanders and Weiser [50] except that the meshes used by the a scheme 
and the latter schemes are all staggered in time. The key features of the a scheme that 
distinguish it from the latter schemes include: (i) the mesh values of both the dependent 
variable and its spatial derivative are considered as the independent variables, to be solved 
for simultaneously ; and (ii) no interpolation or extrapolation techniques are used in the 
construction of the a scheme. Note that the differences between the latter schemes and an 
extension of the a scheme were also clearly spelled out by Huynh [51]. 

This section is concluded with the following remarks: 

(a) The a scheme can be constructed from a different perspective in which both CEs and 
SEs have the shape of a rhombus [2]. In this alternative construction, the differential 
condition Eq. (1.5) is not assumed. Instead it becomes a result of a local flux conserva- 
tion condition and Eq. (1.4). In other words, the a scheme can be constructed entirely 
from flux conservation conditions and the assumption that u (x,t ,J,u.) is linear in x 
and t. 

(b) The a scheme has many non-traditional features. They were discussed in great detail 
in [2]. 

(c) Because there are two independent marching variables at each mesh point € ft, two 
amplification factors appear in the von Neumann stability analysis of the a scheme [2]. 
It happens that these two factors are identical to those of the Leapfrog scheme [52] if 
the latter factors arise from a “correct” von Neumann analysis [2], Note that the main 
Leapfrog scheme (excluding its starting scheme which relates the mesh variables at 
the first two time levels), the Lax scheme [52], and the main DuFort-Frankel scheme 
[52] share one special property, i.e., a solution to any one of these schemes is formed by 
two decoupled solutions. Traditionally the von Neumann analysis for these schemes 
is performed without taking into account this decoupled nature. It is explained in 
[2] why such an erroneous analysis will result in a dispersive property prediction that 
makes the dispersion appear worse than it really is. Moreover, because (i) the a 
scheme and the Leapfrog scheme share the same amplification factors, and (ii) the a 
scheme is a two-level scheme while the Leapfrog scheme is a three-level scheme, the a 
scheme can be considered as a more advanced and compact Leapfrog scheme. 
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The fact that the amplification factors of the a scheme are related to those 
of a celebrated classical scheme is only one among a string of similar unexpected 
coincidences encountered during the development of the CE/SE method. As it turns 
out [2,12,25], the amphfication factors of the Lax, the Crank-Nicolson, and the DuFort- 
Frankel schemes also are related to those of some of the extensions of the a scheme. 
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2. Review of the ID Schemes 


In this section, we shall (i) review and reformulate the ID schemes described in [2], 
and (ii) fill a gap in the derivation of Eq. (4.28) in [2]. Not only does the reformulation 
enable the reader to see more clearly the structural similarity between the ID solvers of 
Eq. (1.1) and their Euler counterparts, it also makes it easier for him to appreciate the 
consistency between the construction of the ID CE/SE solvers and that of the 2D solvers 
to be described in the later sections. 


2-1- The a Scheme 


As the first step, the marching procedure of the a scheme will be cast into a form 
slightly different from that given in [2]. To proceed, let the Courant number v = aAt/Ax . 
Also let 

(2.1) 


«)7 def 




for any (j,n) € ft. Hereafter the superscript symbol “+” is used to denote a normalized 
parameter. Using Eq. (2.1), Eqs. (1.6)— (1.8) imply that 


[(i-^+(i-^i; 


[(1 - u)u - (1 



n- 1/2 
1+1/2 


( 2 . 2 ) 


and 


- 1/2 


[(1 + I ')u - (1 - V 2 )u + ] . = [(1 + v)u + (1 - V 2 ) u t] 


(2.3) 


for all (j, n) £ f!. To simplify notation, in the above and hereafter we adopt a convention 
that can be explained using the expression on the left side of Eq. (2.2) as an example, i.e., 


[(1 - u)u + (1 - P 2 )«+]“ = (1 - n)u~ + (1 - v‘)(ul )” 
Moreover, to streamline the future development, we define 

( s +)j + l/2 = [ U ~ (1 + V ) U » J j+1/2 

( <S -)j_l/2 "" [ U + “ U ) U x J j-1/2 

and 

^, a +t n d — f 1 ( « l”- 1 / 2 _ /• \W-l/2 

( U x )j - 2 [( S +U + l/2 ( 5 -lj-l/2j 

By adding Eqs. (2.2) and (2.3) together, and using the above definitions, one has 


(2.4) 

(2.5) 

(2.6) 


• 7=2 


\ (! - Zy )( S +)"+l/2 2 + + ^)( 5 -)j-l/2 


(j,n) e ft (2.7) 
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Let 1 — v 2 0, i.e., 1 — v ^ 0 and 1 + v ^ 0. Then Eqs. (2.2) and (2.3) can be divided by 
(1 — u) and (1 + u), respectively. By subtracting the resulting equations from each other 
and using Eqs. (2.4)-(2.6), one has 

«)” = « + )j\ (j,n) € SI (2.8) 

Because both and are explicit functions of the marching variables 

at the {n — l/2)th time level, Eqs. (2.7) and (2.8) form the explicit marching procedure for 
the a scheme. Note that these equations can be obtained from the inviscid form of the a-p 
scheme, i.e., Eq. (2.14) in [2]. Also note that the superscript symbol “a” in the parameter 
( U x + )J i s introduced to remind the reader that Eq. (2.8) is valid for the a scheme. 

2.2. The a-e Scheme 

In the a-e scheme [2], CE+(y, n) and CE_(y,n), which are depicted in Figs. 4(d) and 
4(e), respectively, are not considered as conservation elements, i.e., Eq. (1.8) is no longer 
applicable. Instead, one assumes that 

f h* • ds = 0, ( j,n ) G 0 (2-9) 

JS(CE(j,n )) 

where CE(jf,n) is the union of CE + (jf,n) and CE_(j,n) (see Fig. 4(f)). In other words, 
CE(j, n) is a conservation element in the a-e scheme. Again the local conservation condition 
Eq. (2.9) leads to a global conservation condition [2], i.e., the total flux of h* leaving the 
boundary of any space-time region that is the union of any combination of new CEs will 
also vanish. 

It was explained in [2] that Eq. (2.7) follows directly from Eq. (2.9). As a result, the 
former is also valid in the a-e scheme. The a-e scheme is formed by Eq. (2.7) and another 
equation that differs from Eq. (2.8) only in the expression on the right side. To show more 
clearly the similarity of the ID schemes and their 2D versions to be described in the later 
sections, in the following, the counterpart of Eq. (2.8) in the a-e scheme will be rederived 
from a perspective different from that presented in [2]. 

Let (j,n) £ 0. Then ( j ± 1/2, n - 1/2) £ fl. Let 

U j± 1/2 = U j± 1/2 + (At/2 ){. u t)j±Ji2 (2.10) 

Substituting Eqs. (1.5) and (2.1) into Eq. (2.10) and using the definition v — aAt/Ax, one 
has 

«£i/ 2 =[»-2r «+];;*// (2.11) 

Note that, by definition, (j ± 1/2, n) ^ n if (i,n) £ Cl. Thus associated with a 

mesh point ^ Cl. The reader is warned that similar situations may occur in the rest of this 
paper. 
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According to Eq. (2.10), ^± 1/2 can be interpreted as a first-order Taylor’s approxi- 
mation of u at ( j i 1/2, n). Thus 


n _ n 

/ c+xn def U j+l/2 U j-l/2 
\ U x )j 4 


AX 

"T 


n _ 7 t*n 
U j+ 1/2 J-l/2 

AX 


( 2 . 12 ) 


is a central- difference approximation of du/dx at (j, n), normalized by the same factor 
ax/4 that appears in Eq. (2.1). Note that the superscript “c” is used to remind the reader 
of the central- difference nature of the term • In the a-e scheme, Eq. (2.8) is replaced 

by (»J)" = K + )" + 2(K + -»;+)" (2.13) 

where e is a real number. 

At this juncture, note that, at each mesh point (j, n) £ fi, Eqs. (2.7) and (2.8) are 
the results of two conservation conditions given in Eq. (1.8). Because Eq. (2.13) does not 
reduce to Eq. (2.8) except in the special case e = 0, at each mesh point (j, n) £ fi, generally 
the a-e scheme satisfies only the single conservation condition Eq. (2.9) rather than the 
two consevation conditions Eq. (1.8). However, because (u x +)J generally is present on the 
right side of Eq. (2.13), the a-e scheme generally will still be burdened with the cost of 
solving two conservation conditions at each mesh point. The exception occurs only for the 
special case e — 1/2, under which Eq. (2.13) reduces to (u^)j — (u x ^~)j . 

Note that the first part of the expression on the right side of Eq. (2.13), i.e., {u ^)™ , 
emerges from the development of the non-dissipative <x scheme. As a result, it is the non- 
dissipative part. On the other hand, the second part, whose magnitude can be adjusted 
by the parameter e, represents numerical dissipation introduced by the difference between 
the central difference term ( u c +)] and the non-dissipative term «+)? - Thus one may 
heuristically conclude that the numerical dissipation associated with the a-e scheme can 
be increased by increasing the value of 6. It was shown in [2] that this conclusion is indeed 
valid in the stability domain of the a-e scheme, i.e., 

0 < e < 1, and v 2 < 1 (2.14) 


According to Eqs. (2.4)-(2.6), (2.11) and (2.12), both «+)? and «+)? are explicitly 
dependent on v (and therefore explicitly dependent on At). However, ( u — u x Jr )j is not 
dependent on v . As a matter of fact, it can be shown that 


(u 


C+ 

X 


U 


a + 



1 \(u + \ n ~' 12 

2 0 + 1/2 


+ «) 


n- 1/2 
3 - 1/2 



n- 1/2 

j+1/2 


— U 


n- 1/2 
3- 1/2 


) 


(2.15) 


Let ( du x )j be the parameter defined by Eq. (3.2) in [2]. Then it can be shown that 


(u 


C+ 

X 


A 7* 

< + )l = T (^.)? 


(2.16) 
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Note that, in the original development [2], ( du x )J was introduced to break the sym- 
metry of the stencil of the a scheme with respect to space-time inversion. This symmetry 
breaking results in the a-e scheme that was originally defined by the matrix equation 
Eq. (3.6) of [2]. Its two component equations are Eq. (2.7) and 


«)? = 


u 


a+\n 

h 


+ e 



71 — 1 /2 

j+1/2 


+ «) 


n-1/2 

i- 1/2 


1 

2 


(u 


n- 1/2 

j+1/2 


— U 


n — 1/2 
1/2 


) 


(2.17) 


with the latter being equivalent to Eq. (2.13). It should be emphasized that the fact that 
( u t)J — ( u x+) 1 j w hen t — 1/2 , and that therefore the a-e scheme can be considered as a 
central- difference scheme in this special case, was a later accidental discovery . 

2.3. The Euler a Scheme 


For a reason that will soon become obvious to the reader, reformulation of the inviscid 
(p = 0) version of the Navier-Stokes solver described in Section 5 of [2] will precede that 
of the Euler solvers described in Section 4 of [2]. Because the inviscid version is also an 
Euler solver and, like the a scheme, is free of numerical dissipation if it is stable, it will be 
referred to as the Euler a scheme. 

To proceed, consider the Euler equations [2] 

^ + ^ = °, m = !,2,3 (2.18) 

where (i) a m , m ~ 1,2,3, are the independent flow variables to be solved for, and (ii) /m, 
m = 1,2, 3, are known functions [2] of u m , m = 1,2,3. Assuming that the physical solution 
is smooth, Eq. (2.18) is a result of the more fundamental space-time flux conservation laws 

<f h m -ds = 0, m = 1,2,3 (2.19) 

JS(V) 

where h m — ( f m . a m ) , nr = - 1,2,3. 

To proceed, let (i) 


fm t k = df m /du k , m,k = 1,2,3 (2.20) 

and (ii) F + be the 3x3 matrix formed by (At/Ax)/ mi fc, m,k = 1,2,3. Note that, as a 
result of (ii), F + = (At/ ax)F where F is the matrix that appears in Eq. (4.8) in [2]. Also 
let (u m )J be the numerical version of u m at any (j, n) 6 ft. Because f m and f m are 
functions of u m , for any ( j,n ) 6 ft, we can define (/ m )” and (/ m ,fc)!> to be the values of 
f m and respectively, when u m , m = 1,2,3, respectively, assume the values of ( u m V>, 

Tn = 1,2,3. Furthermore, because / m , m = 1,2,3, are homogeneous functions of degree 1 
[53, p. 11] in the variables u m , m = 1,2,3, we have 


3 


(/m)j — y/(/m,<;)j t ( u fc)? 

ife=l 


( 2 . 21 ) 
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Note that Eq. (2.21) is not essential in the development of the ID CE/SE Euler solvers. 
However , in some instances, it is used to recast some equations into more convenient forms. 

For any (x,t) G SE(j,n), u m {x,t), f m (x,t ) and h m (x,t) are approximated by 


and 


; j,») d = Ml + - *;) + («-)?(* - n ( 2 . 22 ) 

f( n {x,t;j,n) = (f Tn )l + (f T nx) 1 j(x — x j) + (fmt)j(t-t ) (2.23) 

h* m (x,t-,j,n) = (fn(x,t;j,n), t i* m (x,t]j,n)) (2.24) 


respectively [2]. Here (i) {u m )] and (tt m ,)£ are the independent marching variables to be 
solved for, and (ii) ( f mx )], (fmt)], and (u m t)? are the functions of (u m )J and {u mx )l , 
m = 1,2,3, defined by Eqs. (4.10), (4.11), and (4.17) in [2]. 

For all (j,n) G f 1, we assume that 


( i ) h^- ds = 0, m = 1,2,3 (2.25) 

JS{CE ± (i,n)) 

Note that Eqs. (2.18), (2.19) and (2.25) are the Euler counterparts of Eqs. (1.1), (1.3) 
and (1.8), respectively. With the aid of Eqs. (2.22)-(2.24), Eq. (2.25) implies that, for all 

(j,n)eU, 


Ml 


t \n-l /2 , AS 

\ U m)j±i/2 ^ a 


\n- 1/2 

U mx )j± 1 j 2 



± 


A t 
AX 



n — 1 /2 
j±l/2 


(fm)l 


± 


(aQ 2 

4ax 


+ (/•»<)? = 


0. 


(2.26) 


Eq. (2.26) is the inviscid version of the Navier-Stokes marching scheme originally given in 
Eq. (5.19) of [2]. 

For each (j,n) G 0, let (i) 


u 


+ 


\n def \n 

,)j — . yUmxJj , 


m = 1,2,3 


(2.27) 


(ii) ut and («+)”, respectively, be the 3 x 1 column matrices formed by (u m )" and («+,.)" , 
m = 1,2,3, and^iii) ( F+)" be the 3 x 3 matrix formed by (At/Ax)(/ m ,fc)”, m,k = 1,2,3. 
Then with the aid of Eqs. (4.10), (4.11) and (4.17) in [2], and Eq. (2.21), one can rewrite 
Eq. (2.26) as a pair of matrix equations, i.e. for any ( j,n ) G fl, 

[(/ - F*)u + (I - (i r+ ) 2 ) = [(/ - F+)n- (I - (F + ) 2 )«t]”; i 1 / 2 ! (2-28) 
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and 


[(I + F + )u- (I - (f+) 2 ) 5+]” = [(/ + F+)u + (j ^ (F+) 2 ) S+];- 1 /’ (2.29) 

where / is the 3x3 identity matrix. 

Note that the flux conservation conditions Eqs. (2.2) and (2.3), and its Euler coun- 
terparts, i.e., Eqs. (2.28) and (2.29) share the same algebraic structure. As a matter of 
fact, the former pair will become the latter pair if the symbols 1, v, u and u+ are replaced 
by I, F+, u and u+ , respectively. As a result, Eqs. (2.28) and (2.29) will be solved by a 
procedure similar to that used earlier to extract Eqs. (2.7) and (2.8) from Eqs. (2.2) and 
(2.3). However, because (i) matrix multiplication is not commutative and (ii) the matrix 
(F + )'j i s a function of (n m )j , m = 1,2,3, while v is a simple constant, as will be shown 
shortly, the algebraic structure of the solution to Eqs. (2.28) and (2.29) is more complex 
than that of Eqs. (2.7) and (2.8). 

To proceed, let (j,n) £ ft and 


and 




<*+)?+$ =' (” - (I + 

(2.30) 

(»*- )”-i/2 =' [« + u - 

(2.31) 

. (2.28) and (2.29) implies that 


[{[(/- F*)s + }^ + [</ + f + )r-]“: l 1 /’} 

(2.32) 

is equivalent to Eq. (4.24) in [2]; and (ii) Eqs. (2.30)-(2.32) are 


Equation (2.32) represents the first part of the solution to Eqs. (2.28) and (2.29). 
To obtain the second part, one must assume the existence of the inverse of the matrix 
[I - (F+) 2 ]? for all (j, n) £ ft. In the following, we shall briefly discuss the significance of 
this assumption. 

Let v and c be the fluid speed and sonic speed, respectively. They are known functions 
of u m , m = 1,2,3 [2]. For each ( j,n ) £ ft, let vj and c™, respectively, denote the values of 
v and c when u m , m = 1,2,3, respectively, assume the values of ( u m )J , m = 1,2,3. Let 

(">)? = =W-<5>. M* *£<•? + «?) (2-33) 


Then, by using (i) the relation F+ = (a t/Ax)F, (ii) the fact that the eigenvalues of the 
matrix F are v — c,v and v + c (see Eq. (4.8) in [2]), and (iii) the fact that the eigenvalues 
°f f{A) are /( Ai), /( A 2 ), /(A 3 ), . . ., /(A n ) if the eigenvalues of a matrix A are Aj, A 2 , A 3 , 
. . ., A n and f(A) is a polynomial of A, one concludes that the eigenvalues of [/ — ( F+ ) 2 ]^ 
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are [1 - ((t'/)") 2 ], l = 1,2,3. Because any square matrix is nonsingular (and therefore its 
inverse exists) if and only if all its eigenvalues are nonzero [54, p.14], one concludes that 
the inverse of [ I - (F + ) 2 ]” exists if and only if 

[MilVl. f = 1,2,3 (2.34) 

In this paper, we shall assume a more restrictive condition than Eq. (2.34), i.e., for all 
(j,n) £ ft, the local Courant number uj < 1. Here 

uj d = max{|(i^i)j|, |(^)7|, |(^3 )j 1} (2.35) 


Note that, because 

(/ - F+)(J + F + ) = (. I + F + )(I -F+) = I- (F+) 2 (2.36) 

the inverse of [I ± (F + )}™ exists if the inverse of [I — (F + ) 2 ]” exists. 

Let (j,n) £ f). Let the marching variables at the (to — l/2)th time level be given. Then 
iPj can be evaluated using Eq. (2.32). Because [/ i is a function of Uj , it follows that 

(S+)“ =' [(/ - F + )"]'‘ [(/ - F+)S - (I - ( 2 - 37 > 

(§-)] = [(/ + F + )"] -1 l(I + F + )u+(I-(F + r)i!i]’' j :; > / 2 2 (2.38) 

and 

« + )" = \(S+ - S-)” < 2 - 39 ) 

can also be evaluated. Note that, in the above and hereafter, the inverse of a matrix A is 
denoted by A” 1 . 

To obtain the second part of the solution to Eqs. (2.28) and (2.29), they are multiplied 
from the left by 

[(/-.F+)?]- 1 and [(/ + F+)”r‘ 

respectively. Let the resulting expressions be subtracted from each other. Then, with the 
aid of Eq. (2.36), one obtains 

(t £)7 = (tS + );> 0»en ( 2 - 4 °) 

Equations (2.32) and (2.40) define the marching procedure of the Euler a scheme. Note 
that the superscript symbol “a” in (u® -1- )^ is intoduced to remind the reader that Eq. (2.40) 
is valid for the Euler a scheme. 

It has been shown by numerical experiments that the Euler a scheme is neutrally stable 
in the interior of the computational domain up to at least a thousand time steps when 
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vj < \ for all (j, n) 6 f l . In these numerical experiments involving a shock-tube problem, 
the computational domain was allowed to grow with time, so that the undisturbed fluid 
state could always be prescribed at the computational boundaries as the exact solution. 
As a matter of fact, by using an analysis similar to that given at the end of Sec. 6 in 
[7], one can show that the linearized form of the Euler a scheme is neutrally stable when 
v™ < 1 for all (j, n) € fi. 

The parameters (S+)" and (S-)J can be evaluated by using Eqs. (2.37) and (2.38) di- 
rectly. This direct evaluation involves inverting two 3x3 matrices which is computationally 
costly. In the following, we shall describe a more efficient approach. 

According to Eqs. (2.37) and (2.38), (S + )J and (SL)” are the solutions to 


(I - F+)?(S + )? = [(/ - F + )S- (/ - (F+) J ) u +]’-■/; (2.41) 

and 

(/ + F+)"(5-)J = [(/ + F>-+ (I - (F+) 2 ) ut}^ (2.42) 

respectively. Note that: (i) each of Eqs. (2.41) and (2.42) represents a system of three 
scalar equations; (ii) because of the reason given in the paragraph preceding Eq. (2.37), 
the coefficients of both systems are known if the marching variables at the ( n — l/2)th 
time level are given, i.e., both systems can be considered as linear; and (iii) because of the 
assumption v™ < 1, each system has a unique solution. As a result of (i)-(iii), both (5+)” 
and (5'-)" can be solved efficiently by using the Gaussian elimination method. 

2.4. The Simplified Euler a scheme 

In implementing the Euler a scheme, two systems of linear equations must be solved 
for each ( j,n ) 6 fl. As a result, the Euler a scheme is locally implicit [1, p.22]. In this 
subsection we shall develop a simplified version that is completely explicit. 

To proceed, the expressions 


[(I-F+)”] and [(/ + F+)"] 1 
in Eqs. (2.37) and (2.38) are approximated by 


\l ~ F + )”" I/2 


1 “I 


0 + 1/2 

respectively. As a result, one has 


and [(/ + F+)”-’/ 2 ] 


(5 + )? « (.-+);;# and (s_)“«(i-_);:>/ 2 2 


n- 1/2 


(2.43) 


where (i+JJ+j 1 ^ 2 and (i*-)^^ 2 are defined in Eqs. (2.30) and (2.31), respectively. Let 

def 1 T/-. \n — 1/2 / \n— l/2l /r . .. 

K )j = 2 [( 5 +)/+i/ 2 " (^-)y_i/2 J (2-44) 
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Then (i) {u“' + )] can be evaluated explicitly, and (ii) as a result of Eqs. (2.39) and (2.43), 
Eq. (2.40) can be approximated by 

(*£)? = K' + )”> ( 2 - 45 ) 

The marching procedure defined by Eqs. (2.32) and (2.45) preferred to as the simplified 
Euler a scheme. Note that the superscript symbol V” in («“'+)" is introduced to remind 
the reader that Eq. (2.45) is valid for the simplified Euler a scheme. 

Generally CE±(j,n), (j,n) G ft, are not conservation elements in the simplified 
scheme. However, because Eq. (2.32) is equivalent to the conservation condition [2] 


/ h* m -ds = 0, (j,n) G ft and m = 1,2,3 (2.46) 

Js(CE(j,n)) 

CE(j,n), (j,n) G ft, are the conservation elements in the simplified scheme. 

Note that by replacing the symbols s_, u“ + , u, u+, 1 and u in Eqs. (2.4)-(2.8) by 
s + , s-, u, ut > I and E + , respectively, these equations will become Eqs. (2.30), (2.31), 
(2.44), (2.32) and (2.45), respectively. In other words, the a scheme and the simplified Euler 
a scheme share the same algebraic structure . 

The simplified Euler a scheme generally is unstable. However, as will be shown shortly, 
this scheme can be extended to become the simplified Euler a-e scheme which does have a 
large stability domain. 

2.5. The Euler a-e Scheme 

The process by which the a-e scheme was constructed from the a scheme will be used 
to construct the Euler a-e scheme from the Euler a scheme. 

In the Euler a-e scheme, the conservation conditions given in Eq. (2.46) are assumed. 
Because Eq. (2.32) is equivalent to Eq. (2.46), the former is also a part of of the Euler a-e 
scheme. The Euler a-e scheme is formed by Eq. (2.32) and another equation that differs 
from Eq. (2.40) only in the expression on the right side. 

To proceed, let (j,n) G ft and 

*&/» = «£# + (“/* < 2 - 47 > 


where (u t )J ±1 V 2 2 is the column matrix formed by (tt ro t)J ± ^ rn — 1,2,3. With the aid of 
Eqs. (4.10) and (4.17) in [2], Eq. (2.47) implies that 


r* / n 

u j± 1/2 


= (u-2F+«+) 


n— 1 /2 
j±l/2 


(2.48) 


Let 


,7'n — u ln , 

t \n def U j + l/2 U J-l/2 

\ u x )j — 4 


(2.49) 


NAS AyTM— 1 998-208843 


19 



Then the Euler a-e scheme is formed by Eq. (2.32) and 

= (^ + ); + 2e(^ + - i 5+)? (2.50) 

where e is a real number. Obviously Eq. (2.50) reduces to (i) Eq. (2.40) when e = 0, and 

(“) )j = (“x + )7 when e = 1/2. Also it has been shown numerically that (i) the Euler 
a-e scheme generally is stable if 


0 < e < 1, and uj < 1 for all (j,n) E ft (2.51) 

and (ii) the numerical dissipation associated with the scheme increases as the value of e 
increases. Note that Eqs. (2.47)-(2.50) are the Euler counterparts of Eqs. (2.10)-(2.13), 
respectively. 

2.6. The Simplified Euler a-e Scheme 

According to Eq. (2.50), excluding the special case e = 1/2, implementation of the 
Euler a-e scheme also requires the evaluation of («“+)" and therefore (see Eqs. (2.37)- 
(2.39)) the solution of Eqs. (2.41) and (2.42). Thus the Euler a-e scheme is locally implicit 
if c ^ 1/2. A totally explicit variant, referred to as the simplified Euler a-e scheme, is 
defined by Eq. (2.32) (or, equivalently, Eq. (2.46)) and 

(<£)? = + 2 e(ul+ - ul'+)J (2.52) 

Obviously the simplified Euler a-e scheme (i) reduces to the simplified Euler a scheme 
when e = 0, and (ii) is identical to the Euler a-e scheme when e = 1/2. 

Note that by replacing the symbols , s_, «“+, u, u+ , u' , u c x + , 1 and v in Eqs. (2.4)- 
(2.7) and (2.11)-(2.13) by 5+, s*_, u x + , u, u', , / and F + , respectively, these 

equations will become Eqs. (2.30), (2.31), (2.44), (2.32), (2.48), (2.49) and (2.52) respec- 
tively. In other words, the a-e scheme and the simplified Euler a-e scheme share the same 
algebraic structure. 

It has been shown numerically that the simplified Euler a-e scheme is stable if 


0.03 < e < 1, and vj <1 for all (j,n) E f) (2.53) 

A comparison between Eqs. (2.51) and (2.53) reveals that the simplified version is only 
slightly less stable than the original version. 

According to Eqs. (2.30), (2.31), (2.44), (2.48) and (2.49), both {u c +)j and (u“'+)^ 

are explicitly dependent on the the matrices 311(1 ( F+ )jI$ ( and therefore 

explicitly dependent on At). However, (u£+ - u x + )J is free from this dependency. Let 
(i) (du mx )J be the parameter defined by Eq. (4.26) in [2], and (ii) (du x )^ be the column 
matrix formed by (du mx )J, m = 1,2,3. Then it can be shown that 






n-l/2 

1 + 1/2 



\ (“+fi 1 /2 1/2 2 ) = (2-54) 
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With the above preliminaries, we are now ready to provide a proof for Eq. (4.28) in 
[2]. Note that the last equation was introduced in [2] simply as a “natural generalization 
of Eq. (3.10) in [2]. 

To proceed, note that Eq. (2.47) is the matrix form of Eq. (4.27) in [2], i.e., uj" 1/2 is 
the column matrix formed by KJjii/a* m = 1,2,3, which were introduced in the latter 
equation. As a result, with the aid of Eqs. (2.27), (2.49) and (2.54), Eq. (2.52) can be 
rewritten as 


W = f (0“+./2 - (<0"- l /2 /“ + ( 2e - !)(<*“•».)” 


(2.55) 


i.e., Eq. (4.28) in [2]. 

Because Eqs. (4.24) in [2] are equivalent to Eq. (2.32), the Euler scheme defined by 
Eqs. (4.24) and (4.28) in [2] is identical to the simplified Euler a-e scheme. 

2.7. The a-e-a-p Scheme and Its Euler Versions 

Consider the a-e scheme defined by Eqs. (2.7) and (2.13). If discontinuities are present 
in a numerical solution, the above scheme is not equipped to suppress numerical wiggles 
that generally appear near these discontinuities. In the following, we shall describe a 
remedy for this deficiency. 

Let . 

(2 ’ 56) 

Then it can be shown that 

« + )" - \ [(»;+)?+ «-)”] (2 - 57) 

i.e., «+)" is the simple average of «+)? and ( )?■ Next, let the function W 0 be 
defined by (i) W o (0,0,a) = 0 and (ii) 


W 0 (x_, x+;a) 


X+|“x_ + |x_|“x + 
|x+|“ + |x_|“ 


(|x+| + |x_|>0) 


(2.58) 


where x + , x_ and a > 0 are real variables. Note that (i) to avoid dividing by zero, in 
practice a small positive number such as 10 -60 is added to the denominator in Eq. (2.58), 
and (ii) W 0 (x_,x+; a), a nonlinear weighted average of x_ and x+, becomes their simple 
average if a = 0 or |x_ | = |x+|. Furthermore, let 


«+)? d = W 0 {« + + )],(u c x ±)]-,a) (2.59) 

Note that the superscript “w” is used to remind the reader of the weighted-average nature 
of the term (tt® + )£. With the aid of the above definitions, a more advanced scheme, 
referred to as the a-e-a-fi scheme, can be defined by Eq. (2.7) and 


(«+)? = K 


a+ ^ n + 2e«+ - <+ )? + (3{u™ + 


h 


U 


c+\n 
n 


(2.60) 
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Here (3 > 0 is another adjustable constant. Note that Eq. (2.60) can be rewritten as 


«)] = P W o + (1 -/?)«+)? + (2e - 1 )«+ - <+)? (2.61) 


It can be shown easily that the a-e-a-f3 scheme reduces to the a-e scheme if a = 0 or (3 = 0. 

The expression on the right side of Eq. (2.60) contains three parts. The first part is 
a non-dissipative term («“+)?. The second part is the product of 2e and the difference 
between the central difference term (u c x + )j and the non-dissipative term (u“+)? . The 
third part is the product of (3 and the difference between a weighted average of (u^)™ 
and (u c x t)” and their simple average. Numerical dissipation introduced by the second part 
generally is effective in damping out numerical instabilities that arise from the smooth 
region of a solution. But it is less effective in suppressing numerical wiggles that often occur 
near a discontinuity. On the other hand, numerical dissipation introduced by the third 
part is very effective in suppressing numerical wiggles. Moreover, because the condition 
\(u c x \)‘j | = |(u*l)"| more or less prevails and thus the weighted average is nearly equal 
to the simple average (see the comment given immediately following Eq. (2.58)) in the 
smooth region of the the solution, numerical dissipation introduced by the third part has 
very slight effect in the smooth region. 

From the above disscusion, one concludes that there are two different types of numer- 
ical dissipation associated with the a-e-a-/3 scheme and they complement each other. As 
a result, the a-e-oc-f3 scheme can handle both small disturbances and sharp discontinuies 
simultaneously if the values of e, a and f3 are chosen properly (usually e = 1/2, a = 1,2 
and f3 = 1). Also note that, to give the CE/SE method more flexibility in controlling local 
numerical dissipation, the parameters e and /? can even be considered as functions of local 
dynamical variables and mesh parameters (see Sec. 8). 

Similarly, the Euler a-e scheme and the simplified Euler a-e scheme can be modified 
to become the Euler a-e-ai-f3 scheme and the simplified Euler a-e-oc-f3 scheme, respectively, 
by simply replacing Eqs. (2.50) and (2.52) with 


«)" = K + )" + 2e(ul+ - «;+)"+/?«+ - ^ + )” 


and 


(2.62) 

(2.63) 


«)" = «' + )” + 2e(«T - K ,+ )7 +£« + - )” 

respectively. Here (u ™ +)” is the 3x1 column matrix formed by 

rn= 1,2,3 

where 

« + ,±)" ±l(«)" ± , /2 ~ (« m )") (2.64) 

with an d (u m )™ being the mth components of j^ 2 and respectively. 
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2.8. The ID CE/SE Shock-Capturing Scheme 

Let e = 1/2 and (3 = 1. Then the Euler a-e-a-(3 scheme and the simplified Euler 
a-t-a-fd scheme reduce to the same scheme. The reduced scheme is defined by Eq. (2.32) 
and 

= m=1 > 2 * 3 (2 ' 65) 

where (j,n) (E fl. 

The above scheme is one of the simplest among the Euler solvers known to the authors. 
The value of a is the only adustable parameter allowed in this scheme. Because it is totally 
explicit and has the simplest stencil, the scheme is also highly compatible with parallel 
computing. Furthermore, it has been shown that the scheme can accurately capture shocks 
and contact discontinuities with high resolution and no numerical oscillations. For these 
distinctive features and for convenience of future reference, the reduced scheme will be 
given a special name, i.e., the ID CE/SE shock-capturing scheme. Note that this scheme 
with a = 1 is implemented in the two shock-tube solvers referred to in Sec. 1. Consider 
only the case that all spatial boundary points (j, n) E fi are at the time levels n = 0, 1, 2, . . . 
(see Fig. 4(a)). The non-reflecting boundary conditions used in the first solver, i.e., the 
one listed in Appendix A, are: (i) 


sf = \% and (> = "=1,2,3,... (2.66) 

if (j,n) is a mesh point on the right spatial boundary; and (ii) 

u] = *j+l/2 and )j = (“£ )"+ 1/2 ’ n = 1, 2, 3, . . . (2.67) 

if (j,n) is a mesh point on the left spatial boundary. On the other hand, for the alternate 
solver, the steady-state boundary conditions 


Uj = Uj and (^)j=(^)ji n- 1,2,3,... 
is imposed at any mesh point (j, n) on the right or left spatial boundary. 


( 2 . 68 ) 
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3. Geometrical Description of Conservation Elements 
in Two Spatial Dimensions 

In Sec. 2, it was established that, for each ID CE/SE solver, there were 2 M indepen- 
dent marching variables per mesh point with M being the number of conservation laws to 
be solved. Because M conservation conditions are imposed over each CE, two CEs were 
introduced at each mesh point such that both the ID a scheme and the ID Euler a scheme 
can be constructed by solving, at each mesh point ( j,n ) € 0, for the 2 M variables using 
the 2 M conservation conditions imposed over CE_(j, n) and CE+(y,n). 

As will be shown in the following sections, for each 2D CE/SE solver, there are 3 M 
independent marching variables per mesh point. As a result, construction of the 2D a 
scheme and the 2D Euler a scheme demands that three CEs be defined at each mesh 
point. In this section, only the basic geometric structures of these CEs will be described. 

Consider a spatial domain formed by congruent triangles (see Fig. 5). The center 
of each triangle is marked by either a hollow circle or a solid circle. The distribution of 
these hollow and solid circles is such that if the center of a triangle is marked by a solid 
(hollow) circle, then the centers of the three neighboring triangles with which the first 
triangle shares a side are marked by hollow (solid) circles. As an example, point G , the 
center of the triangle A BDF, is marked by a solid circle while points A, C and E , the 
centers of the triangles A FMB } A BJD and A DLF, respectively, are marked by hollow 
circles. These centers are the spatial projections of the space-time mesh points used in the 
2D CE/SE solvers. 

To specify the exact locations of the mesh points in space-time, one must also specify 
their temporal coordinates. In the 2D CE/SE development, again we assume that the 
mesh points are located at the time levels n = 0, ±1/2, ±1, ±3/2, . . . with t = nAt at the 
nth time level. Furthermore, we assume that the spatial projections of the mesh points at 
a whole-integer (half-integer) time level are the points marked by hollow (solid) circles in 
Fig. 5. 

Let the triangles depicted in Fig. 5 lie on the time level n = 0. Then those points 
marked by hollow circles are the mesh points at this time level. On the other hand, those 
points marked by solid circles are not the mesh points at the time level n = 0. They are 
the spatial projections of the mesh points at half-integer time levels. 

Points A, C and E , which are depicted in Figs. 5 and 6(a), are three mesh points at 
the time level n = 0. Point G", which is depicted in Fig. 6(a), is a mesh point at the time 
level n = 1/2. Its spatial projection at the time level n = 0 is point G. Because point G 
is not a mesh point, it is not marked by a circle in the space-time plots given in Figs. 6(a) 
and 6(c). Hereafter, only a mesh point, e.g., point G' , will be marked by a solid or hollow 
circle in a space-time plot. 

The conservation elements associated with point G' are defined to be the space-time 
quadrilateral cylinders GF ABG' F' A! B ' , GBCDG' B'C'D', and GDEFG'D'E'F' that are 
depicted in Fig. 6(a). Here (i) points B , D and F are the vertices of the triangle with 
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point G as its center (centroid) (see also Fig. 5), and (ii) points A ' , B', C', D ', E' and F 
are on the time level n = 1/2 with their spatial projections on the time level n 0 being 
points A , B, C, D, E and F, respectively. 

Point G' is a mesh point at a half-integer time level. For a mesh point at a whole- 
integer time-level, the conservation elements associated with it can be constructed in a 
similar fashion. As an example, consider Fig. 6(b). Here points B 1 (B"), /'(/"), 7' (J"), 
K' ( K "), D' (!>"), G' (G") and C' (C") are on the time level n = 1/2 (n = 1) with their 
spatial projections on the time level n = 0, respectively, being the points B , /, J, K , D , 
G and C that are depicted in Fig. 5. Point C" is a mesh point at the time level n = 1. By 
definition, the conservation elements associated with point C" are the quadrilateral cylin- 
ders C'J'K , D'C"J"K"D", C'D , G'B'C"D"G"B" and C' B ' /' J'C" B" I" J" . The relative 
space-time positions of the six CEs associated with mesh points G' and C" are depicted 
in Fig. 6(c). 

Recall that, in the development of the ID a scheme, a pair of diagonally opposite 
vertices of each CE±(j,n) (see Figs. 4(d) and 4(e)) are assigned as mesh points. Fur- 
thermore, the boundary of each CE-t(j, n) is a subset of the union of the SEs associated 
with the two diagonally opposite mesh points of this CE. In the 2D development, as seen 
from Figs. 6(a)-(c), two diagonally opposite vertices of each CE are also assigned as mesh 
points. In Sec. 4, we shall define the SEs such that even in the 2D case, the boundary of a 
CE is again a subset of the union of the SEs associated with the two diagonally opposite 
mesh points of this CE. 

As a preliminary to the derivation of several equations to be given in Sec. 4, this 
section is concluded with a discussion of several geometric relations involving point G and 
the vertices of the hexagon ABCDEF that are depicted in Fig. 5. By using the facts that 
(i) points A, C, E and G are the geometric centers of four neighboring congruent triangles 
A FMB, A BJD, A DLF and A BDF, respectively; and (ii) any two of the above four 
triangles form a parallelogram (note: two congruent triangles sharing one side may not 
form a parallelogram), it can be shown that: 

(a) CD, CF , DG and AF are parallel line segments of equal length. 

(b) AB, GC, FC and ED are parallel line segments of equal length. 

(c) DC, GD , AG and FE are parallel line segments of equal length. 

(d) Point G is the geometric center of the hexagon ABC DE F and the triangle AC E. 

Note that the line segments GA, CC,GEAC,CE and EA are not shown in Fig. 5. Also 
note that, because the hexagon BIJKDG (depicted in Fig. 5) is congruent to the hexagon 
ABCDEF , a set of geometric relations similar to those listed above also exists for the 
vertices and the center of the hexagon BIJKDG. 
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4. The 2D a Scheme 


In this section, we consider a dimensionless form of the 2D convection equation, i.e., 


du du du 

dt +a ‘fc +a, dy = ° 


(4.1) 


where a x , and a y are constants. Let x^ — x, x 3 = y, and X 3 — t be the coordinates of a 
three-dimensional Euclidean space E 3 . By using Gauss’ divergence theorem in the space- 
time E 3 , it can be shown that Eq. (4.1) is the differential form of the integral conservation 
law 

<f h ■ ds = 0 (4.2) 

JS(V ) 

Here (i) S'(V') is the boundary of an arbitrary space- time region V in E 3 , (ii) 


r def 
n = 


( a x u , a y u, u) 


(4.3) 


is a current, density vector in E 3 , and (iii) ds — den with dc and n, respectively, being the 
area and the outward unit normal of a surface element on ^(F). It was shown in Sec. 3, 
that E 3 can be divided into nonoverlapping space-time regions referred to as conservation 
elements (CEs). 

In the following analysis, the nontraditional space-time mesh that was sketched in 
Sec. 3 will be rigorously defined. To proceed, the spatial projections of the mesh points 
depicted in Fig. 5 are reproduced in Fig. 7. Note that the dashed lines that appear in 
Fig. 7 are the spatial projections of the vertical interfaces (see Fig. 6(a)-(c)) that separate 
different CEs. Also note that, as a result of the geometric relations listed at the end of 
Sec. 3, any dashed line can point only in one of three different fixed directions. We assume 
that the congruent triangles depicted in Fig. 5 are aligned such that one of the above fixed 
directions is the ar-direction. 

Each mesh point marked by a solid or hollow circle is assigned a pair of spatial indices 
0*> according to the location of its spatial projection. Obviously, a mesh point can 
be uniquely identified by its spatial indices (j, &) and the time level n where it resides. 
According to Figs. 8 and 9, the spatial projections of the mesh points that share the same 
value of j (&) he on a straight line on the x-y plane with this straight line pointing in the 
direction of the k- ( j -) mesh axis. 

Let 

def 

t n = riAi, n = 0, ±1/2, ±1, ±3/2,... (4.4) 

Let j and k be spatial mesh indices with j,k = 0, ±1/3, ±2/3, ±1,.... Let flj denote 
the set of mesh points (j,k,n) with j,k = 0, ±1, ±2, — , and n = ±1/2, ±3/2, ±5/2,.... 
These mesh points are marked by solid circles. Let O 2 denote the set of mesh points 
(j,k, n ) with j,k = 1/3, 1/3 ± 1 , 1/3 ± 2, . . ., and n = 0, ±1, ±2, — These mesh points 
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are marked by hollow circles. The union of fli and fl 2 will be denoted by fl. Note that 
the same symbol fl was also used to denote the set of mesh points used in the ID solvers 
(see Sec. 2). Hereafter, unless specified otherwise, the new definition of fl is assumed. 

Each mesh point (j, k,n) £ fl is associated with (i) three conservation elements (CEs), 
denoted by CE r (j,fc,»), r = 1,2,3 (see Figs. 10(a) and 11(a)); and (ii) a solution element 
(SE), denoted by SE(j, k, n ) (see Figs. 10(b) and 11(b)). Each CE is a quadrilateral cylinder 
in space-time while each SE is the union of three vertical planes, a horizontal plane, and 
their immediate neighborhoods. Note that the CEs and the SE associated with a mesh 
point (j, k,n) £ fli differ from those associated with a mesh point ( j,k,n ) £ fl 2 in their 
space-time orientations. 

By using the geometric relations listed at the end of Sec. 3, one can conclude that 
the spatial coordinates of the vertices of the hexagon ABCDEF , which appears in both 
Figs. 10(a) and 11(a), are uniquely determined by three positive parameters w, b and h 
(see Fig. 12(a)) if (i) one assumes that DA is aligned with the x-direction, and (ii) the 
spatial coordinates of point G (the centroid of the hexagon) a re gi ven. Note th at w, b and 
h , respectively, are the lengths of the line segments DM, MH and BH with (i) DM being 
a median of the triangle A BDF, and (ii) points G, M and H being on the line segment 
DA. Also note that a dashed line in Fig. 7 may appear in other figures as a solid line. 

According to Fig. 7, E 3 can be filled with the CEs defined above. Moreover, it is 
seen from Figs. 10(a), 10(b), 11(a), and 11(b) that the boundary of a CE is formed by the 
subsets of two neighboring SEs. 

Let the space-time mesh be uniform, i.e., the parameters At, w, b, and h are constants. 
Let x j>k and y Jifc be the x- and y- coordinates of any mesh points ( j,k,n ) £ fl. Let x 0 ,o = 0 
and yo,o = 0. Then information provided by Figs. 12(a) and 12(b) implies that 

Xj,k = ( j + k)w + (k - j)b, y jt k = {k - j)h (4.5) 


Let n l5 n 2 , n 3 , n 4 , n 5 , and n 6 be the vectors depicted in Fig. 12(a). They he on the x-y 
plane and are the outward unit normals to AB, BC , CD, DE, EF, and FA, respectively. 
It can be shown that 


(h,-b + w/ 3,0) 
\A 2 + (6 - it>/3) 2 ’ 


U4 = —tli 


n 2 = (0,1,0), n 5 = -n 2 


(4.6a) 

(4.66) 


and 


(— h, b -(- iu/3, 0) 
y/h*+{b + w/3)*' 


h 6 = — n 3 


(4.6c) 


For any (x,y,t) £ SE (j,k,n), u(x,y,t) and h(x,y,t), respectively, are approximated 


by 


U 


'{x,y,t;j,k,n) d =u 7 j ik +{u x )'j' k (x — Xj t k) + (uy)‘j yk {y Vj,k) + ( u t)j,k(^ * ) 
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and 


h*{x,y,t-,j,k,n) = [, a x u*(x,y,t‘J,k,n),a y u*{x,y,t-,j,k,n),u*(x,y,i-,j,k,n )] (4.8) 


where uj k , (“1)^1 (^s i)yk^ an< ^ ( u t)yk are constants within SE(j, k,n). The last four 
coefficients, respectively, can be considered as the numerical analogues of the values of 
u, du/dx , du/dy, and du/dt at (*j,fc,J0',fc,f n ). As a result, the expression on the right 
side of Eq. (4.7) can be considered as the first-order Taylor’s expansion of u(x,y,t) at 
( x j,k,yj,k,t n )• Also note that Eq. (4.8) is the numerical analogue of Eq. (4.3). 

We shall require that u = u*(x,y,t;j,k,n ) satisfy Eq. (4.1) within SE(j,*,n). As a 
result, 

( u t)j,k = ~ [ a x{ux)yk + a y( u y)J,k] (4-9) 

Substituting Eq. (4.9) into Eq. (4.7), one has 


u { x iVi k,n) — Uj k + (u x )y k [( x — x i,k) ~ Oi(f — f n )] 

( u y)j,k [(?/ — Vj,k) — a.y(t — t )] . 


(4.10) 


Thus there are three independent marching variables, i.e., u™ k , (u x )^ k , and {u y ) r > k as- 
sociated with a mesh point ( j,k,n ) G f2. For any ( j,k,n ) G fti, these variables will be 
determined in terms of those associated with the mesh points ( j + 1/3, k + 1/3, n - 1/2), 
(j - 2/3, k + 1/3, n - 1/2), and (; + 1/3,* - 2/3, n - 1/2) (see Fig. 13(a)) by using the 
three flux conservation relations 


£ 


S(CE r (j,k,n)) 


h* ■ ds = 0, r — 1,2,3 


(4.11) 


Similarly, the marching variables at any (j, k,n) G H 2 are determined in terms of those 
associated with the mesh points {j - 1/3, * - 1/3, n - 1/2), ( j + 2/3, * - 1/3, n - 1 /2), and 
(j ~ 1/3,* + 2/3, n — 1/2) (see Fig. 13(b)) by using the three flux conservation relations 
Eq. (4.11). Obviously, Eq. (4.11) is the numerical analogue of Eq. (4.2). 

As a result of Eq. (4.11), the total flux leaving the boundary of any CE is zero. 
Because the flux at any interface separating two neighboring CEs is calculated using the 
information from a single SE, the flux entering one of these CEs is equal to that leaving 
another. It follows that the local conservation conditions Eq. (4.11) will lead to a global 
conservation condition, i.e., the total flux leaving the boundary of any space-time region 
that is the union of any combination of CEs will also vanish. 

In the following, several preliminaries will be given prior to the evaluation of Eq. (4.11). 
To proceed, note that a mesh line with j and n being constant or a mesh line with * and n 
being constant is not aligned with the aj-axis or the y- axis. We shall introduce a new spatial 
coordinate system (Cy) with its axes aligned with the above mesh lines (see Fig. 12(c)). 


NASA/TM — 1 998-208843 


28 



Let e x and e y be the unit vectors in the x- and the y- directions, respectively. Let 

e ^ and e v be the unit vectors in the directions of DF and DB (i.e., the j- and the k- 
directions-see Figs. 12 (a)-(c)), respectively. It can be shown that 

e<; = [(u> - b)e x - he y ] / a( (4-12) 


and 


e, = [("’ + b)e x + he y ]/ AT} 


(4.13) 


where 


and 


d = \DF\ = \J {v> ~ b ) 2 + h 2 

AT) d = \DB\ = y/{w + 6) 2 + h 2 


(4.14) 

(4.15) 


Let the origin of (x,y) also be that of Then, at any point in E 3 , the coordinates 

((, 77 ) are defined in terms of (x,y) using the relation 


(e<; + 1]^ = xe x + ye y 


(4.16) 


Substituting Eqs. (4.12) and (4.13) into Eq. (4.16), one has 


and 


Here 


and 




m itcf 


/ w-b 

h 


w + b \ 
AT) 

h 

&T] / 



(w + b)A( 

2u> 

2wh 

AT] 

V 2w 

(w — b)Ai] 

2wh 


Note that the existence of T -1 , the inverse of T, is assured if wh ^ 0. 


(4.17) 


(4.18) 


(4.19) 


(4.20) 
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With the aid of Eqs. (4.5), (4.18), and (4.20), it can be shown that the coordinates 
((, 77 ) of any mesh point (j,fc,n) 6 Q are given by 


C = j a£, and Tj — k Arj 


(4.21) 


i.e., A< and Ai] are the mesh intervals in the (- and the 77- directions, respectively. 

Next we shall introduce several coefficients that are tied to the coordinate system 
(<,77). Let 


/ a C \ def rp - 1 
\ ) 



(4.22) 


Also, for any ( j,k,n ) € fl, let 


(Ml* 


/ \ 
\( u y)j,fc / 


(4.23) 


where T* is the transpose of T. For those who are familiar with tensor analysis [55], the 
following comments will clarify the meaning of the above definitions: 

(a) (a^,a v ) are the contravariant components with respect to the coordinates (<£, 77 ) for 
the spatial vector whose x- and y- components are a x and a y , respectively. 

(b) ((«c)j,fc’(' u '?)j,fc) are the covariant components with respect to the coordinates ((, 77 ) 
for the spatial vector whose x- and y- components are (u x )J k and (u y )J fc , respectively. 

(c) Because the contraction of the contravariant components of a vector and the covariant 
components of another is a scalar, Eq. (4.9) can be rewritten as 


( u t)j,k — + a v( U v) 1 j,k] (4.24) 

(d) Under the linear coordinate transformation defined by Eqs. (4.17) and (4.18), (£ — 
j a£,77 — kArj) are the contravariant components with respect to the coordinates (£, 77 ) 
for the spatial vector whose x- and y- components are x — xj^ and y—yj,k, respectively. 
Using the same reason given in (c), Eq. (4.10) implies that 


u*(x,y,t;j,k,n) — j,k,n) 


(4.25) 


where 

77*(C,T7,t;j,fc,77) ¥ 7 ll k + [(( - jAt) - a ( (f - *")] 

+ [(>) “ *•“?) - a,(l - <“)] 


(4.26) 


Note that Eqs. (4.24) and (4.25) can also be verified directly using Eqs. (4.18), 
(4.22), and (4.23). 

Next, let (i) 


def 3 a t 

V < = 2AC a ° 


and 




def 3A/ 
2A77 


(4.20), 


(4.27) 
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and (ii) 


(4.28) 


«)i.i 


def 


A( 


Ml 


k 5 


and 


/ -f \n dff ^^7 / \n 

( U i? )j\fc — c K)^ 


The coefficients defined in Eqs. (4.27) and (4.28) can be considered as the normalized 
counterparts of those defined in Eqs. (4.22) and (4.23). Furthermore, because A< and at? 
are the mesh intervals in the £- and rj- directions, respectively, Eq. (4.27) implies that 
(2/3)u c and (2/3 )u v , respectively, are equal to the Courant numbers in the <- and 77 - 
directions, respectively. 

Furthermore, let 


(1)± def n 

<r\i = l-v<-u n 

(4.29) 

o\\ ]± d = ±(1 -u c - u v )( 1 + u ( ) 

(4.30) 

eSS* = ±(1 - «-< - ■'„)(! + 

(4.31) 

(1)± def 1 
(T 21 — 1 + ^ 

(4.32) 

di )±d = ^(l + ^K 2 -^) 

(4.33) 

^3^ d = ±(1 + f'cXl + V V ) 

(4.34) 

(1)± def 1 . 

*31 = 1 + u n 

(4.35) 

^(IJidef ±(1 + ^)(1 + ^) 

(4.36) 


(4.37) 

<rJ5 )± d = 1 + i/c + v v 

(4.38) 

a (2)± def + ^ + ^)(1 - V< ) 

(4.39) 

0*3* d = =F(1 + *>c + * / v )( 1 - 

(4.40) 

= f 1 - 

(4.41) 

oQ* * ±(1 - v ( ){2 + v c ) 

(4.42) 

viV* d = =F(1 - ^c)(! - 

(4.43) 

(2)± def 1 

<t 31 — 1 JA) 

(4.44) 

4?* =' Td - - n) 

(4.45) 

oi-. 11 =' ±(1 - *'.,)( 2 + "») 

(4.46) 
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Note that: 

(a) Each of Eqs. (4.29)-(4.46) represents two equations. One corresponds to the upper 
signs while the other, to the lower signs. 

(b) The definitions given in Eqs. (4.29)-(4.37) will be used in the first marching step of 
the 2D a scheme; while those given in Eqs. (4.38)-(4.46) will be used in the second 
marching step. It is seen that the expressions on the right sides of the former can 
be converted to those of the latter, respectively, by reversing the and ” signs. 
Moreover, for every pair of r and s (r,s = 1,2,3), o-^ )_ and <rlV~ are converted to 
°ra and (Trs* , respectively, if and v-q are replaced by — and — u v , respectively. 


(c) We have 

(«)± , (s)± , _(9)± _ o n 1 9 

°hi ' ° 2 i i ° 3 i — d, q — 1,2 

(4.47) 

and 

cr (?)± i _(9)± i _(?)± 
a 12 ' <t 22 < cr 32 

— -I- 4 - _ n n i n 

— <t 13 + <T 23 + ^33 — 0 , 9 = 1,2 

(4.48) 


To simplify the following development, let 


(j,fc;l,l) d = j + 1/3,* + 1/3 

(4.49a) 

(i, ^ 5 1, 2) = f j — 2/3, k -)- 1/3 

(4.496) 

(j,M,3)= f j + l/3,fc-2/3 

(4.49c) 

(j,fc;2,l) d 4 f j - l/3,*s — 1/3 

(4.50a) 

k] 2, 2) d = j + 2/3, k — 1/3 

(4.506) 

(j,*;2,3) d =?: f j - 1/3, lb + 2/3 

(4.50c) 


Note that (i) (j,k;l,r), r — 1,2,3, are the spatial mesh indices of points A, C , and E 
depicted in Fig. 10(a), respectively, (ii) (j,fc;2,r), r = 1,2,3, are the spatial mesh indices 
of points D, F, and B depicted in Fig. 11(a), respectively, and (iii) the mesh indices on the 
right sides of Eqs. (4.49a,b,c) can be converted to those in Eqs. (4.50a,b,c), respectively, 
by reversing the “+” and ” signs. 

Equation (4.11) is evaluated in Appendix B. Let (j,k,n) £ Q, q with q = 1,2. Then, 
for any r = 1,2,3, the result of evaluation can be expressed as: 


.(«)+ 


<«) + „+ _L_ ~(9) + „,+ 


u + <2 u i + a 


r3 


U ' 


j,k 


r<?- 

rl 


u + a 


(<?)- 

r2 


u~^ + a[ q ^~ 


r3 


U 


+ 


n- 1/2 
(j*k;q,r) 


(4.51) 
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According to Eqs. (4.29)-(4.31), and contain a common factor 

(1 _ U( . — Uy,). Similarly, each of three consecutive pairs of coefficients defined in Eqs. (4.32)- 
(4.46) also contain a common factor. As a result, if one assumes that (i) 1 - - u v ± 0, 

(ii) 1 + u c ± 0, (iii) 1 + ^/0, (iv) 1 + i/< + v v ± 0, (v) 1 - u c ^ 0 and (vi) 1 - */„ / 0, i.e., 

[1 - (i/< + u v ) 2 } (1 - v]) (1 - v*) / 0 (4.52) 

then the six equations (q = 1,2 and r = 1,2,3) given in Eq. (4.51) can be simplified as 



[u + {1 + vc)u+ + (1 + v v )u+Y^ k = 

S (1) 
*1 ’ 

(j,fe,n) G Sll 

(4.53) 


[« - (2 - v<;) u + 4 (1 + k = 

S (1) 

(j,A;,n) G fli 

(4.54) 


|u + (1 4 — (2 — I/ 7j) w ^] . k = 

s (1) 

*3 ’ 

(y,fe,n) G fli 

(4.55) 


[u - (1 - *'<)«< - U “ u v) u n\ j k = 

■S (2) 
*1 J 

{j,k,n) G n 2 

(4.56) 


^■u + (2 4 ^)u^“ — (1 — ^ = 

s (2) 

(j,fc,n) G f^2 

(4.57) 

and 

-in 

« — (1 - ^c) u c + + u v) u n . . = 

L S 1 3,k 

-S (2) 

(j,k,n) G fl 2 

(4.58) 

respectively. 

Here 





Sj 1 * *= f U — (1 4 V()u^ — (1 4 ^jj)^] 

n — 1 /2 

(j,fe,n) G 12i 

(4.59) 


4^ *= f u 4 (2 — — (1 4 Vri)V-^ 

n — 1/2 

(j, fc,n) G fii 

(4.60) 


5 ( 3 a) d = u - (1 4- u ( )u^ + (2 - v v )u+ 

n — 1/2 
(j,fc;l,3) ’ 

(j,fc,n) G fii 

(4.61) 


s[ 2) d = [u 4(1- »d u t + (! “ 

n — 1/2 

(j,fc;2,i) ’ 

(j, fc,n) G 0 2 

(4.62) 


^(2) u — (2 4 I^)u/ 4 (1 — v v) U v 

I n — 1/2 
i(j.fc;2,2) ’ 

(j,k,n) G f *2 

(4.63) 

and 

4^ ^ U 4 (1 — U () U ( - (2 4 V V ) U V 

1 n-1/2 

J(j,fc;2,3) ’ 

{j,k,n) G 0 2 

(4.64) 
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The current 2D a scheme will be constructed using Eqs. (4.53)-(4.58) without assuming 
Eq. (4.52). Note that Eqs. (4.53)-(4.58) imply Eq. (4.51) for any and u v . However, the 
reverse is false unless Eq. (4.52) is assumed. 

Note that the expressions within the brackets in Eqs. (4.53)-(4.55) and (4.59)-(4.61), 
respectively, can be converted to those in Eqs. (4.56)-(4.58) and (4.62)-(4.64) by reversing 
the and ” signs. 

It can be shown that Eqs. (4.53)— (4.55) are equivalent to 


•s* = i 

(1 - l/£ — U v )s[^ + (1 + + (1 + ^)^3^ 

(4.65) 


(«<)’.* = = \ (4 1 * - 4 11 ) 

(4.66) 

and 

«)",» = « + )h = i (4 11 - 4 11 ) 

(4.67) 

where (j, k,n) £ ft*. Also Eqs. (4.56)-(4.58) are equivalent to 


= l 

j,k 3 

(1 + + V v )s { ^ + (1 - ^)^2 2) + (1 - I / t 7 )53 2) J 

(4.68) 


(«<)’.* = K + )".* = § (4 21 - 4 2> ) 

(4.69) 

and 

«)".* = K + )“.» = §(4 2, -4 2) ) 

(4.70) 

where ( j,k,n ) £ fl 2 - 




At this juncture, it should be emphasized that Eqs. (4.65) and (4.68) can be derived 
directly from Eq. (4.51). As a matter of fact, with the aid of Eqs. (4.47) and (4.48), we 
can obtain Eq. (4.65) (Eq. (4.68)) by summing over the three equations with q = 1 (g = 2) 
and r = 1,2,3 in Eq. (4.51). 

The 2D a scheme is formed by repeatedly applying the two marching steps defined by 
Eqs. (4.65)-(4.67) and Eqs. (4.68)— (4.70), respectively. It has been shown numerically that 
it is of second order in accuracy for u? k , (u c )? fc and (u v )] ik assuming that v ( and u v are 
held constant in the process of mesh refinement (note: as a result of Eq. (4.28), the 2D a 
scheme is third order accurate for («+)”*. and (u+)£ fc ). Note that the superscript symbol 
“a” in (u£ + )? fc and (u^ + )£ fc is introduced to remind the reader that Eqs. (4.66), (4.67), 
(4.69) and (4.70) are valid for the 2D a scheme. Although the 2D a scheme is constructed 
using a procedure very much parallel to that used to construct the ID a scheme, the former 
is more complex than the latter in many aspects. One key difference between these two 
schemes is that the 2D a scheme is formed by two distinctly different marching steps while 
the ID a scheme is formed by repeatedly applying the same marching step defined by the 


N ASA/TM — 1 998-208843 


34 



inviscid version of Eq. (2.14) in [2], It is this difference that, in the 2D case, makes it 
necessary to divide the mesh points into two sets fij and D 2 . 


As a preliminary for the stability analysis of the 2D a scheme given in Sec. 6, for any 
(j, fc,n) G D, let 


d = 


u 


+ 


(4.71) 


< ' j,k 


Furthermore, let the six 3x3 matrices Qr q \ q = 1?2, and r = 1,2,3, respectively, be 
the special cases of those defined in Eqs. (5.18)— (5.23) (see Sec. 5) with 6 — 0. Then 
Eqs. (4.65)-(4.70) can be expressed as 

3 

q(j,k,n) = ^Qi q) q[{j,k;q, r ),n - 1/2), (j,k,n) G (4.72) 

r— 1 

Combining Eqs. (4.72) and (4.49a)-(4.50c), one has (i) 

q(j,k,n) = Q^Q^qU + 1 ,M - 1) + Q^Q^q (j, k + 1 ,n - 1) 

+ Q^Q^qU - 1 - 1) + Q^Q^qti - l ,* + - 1) (4<73) 

+ _ ~ i) + + - 1 ,n - 1) 

+ (Q[ 1) Q? ) + q[ 1] Q[ 2) + Qi 1) Qi 2) )q(j,k, n - 1 ) 

where (j,fc,n) G Hi; and (ii) 

q(j,k,n) = (j - l,A:,n - 1) + Q^Qi^qU^- l,n - 1) 

+ + l,A:,n- 1) + Q ( 2 2) Q ( 3 1) 9(i + l,fc - l,n- 1) 

+ Qi 2) Q\ 1) q(j,k + 1 ,» - 1) + Ql^Q^qU - M + 1,« - 1) 

+ (QWi 0 + + Qi 2) Q ( s 1) )9(i,fc,« - 1) 


where (j,fc,n) E 0 2 - Note that (i) Eq. (4.73) relates the marching variables at two adjacent 
half-integer time levels; and (ii) Eq. (4.74) relates the marching variables at two adjacent 
whole-integer time levels. 

The 2D a scheme has several nontraditional features. They are summarized in the 
following comments: 

(a) As in the case of the ID a scheme, the 2D a scheme also has the simplest stencil 
possible, i.e., in each of their two marching steps, the stencil is a tetrahedron in 3D 
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space-time with one vertex at the upper time level and the other three vertices at the 
lower time level. 

(b) As in the case of the ID a scheme, each of the six flux conservation conditions asso- 
ciated with the 2D a scheme., i.e., those given in Eq. (4.51) (q = 1,2 and r = 1,2,3), 
represents a relation among the marching variables associated with only two neigh- 
boring SEs. 

(c) As in the case of the ID a scheme, the 2D a scheme also is non-dissipative if it is 
stable. It is shown in Sec. 7 that the 2D a scheme is neutrally stable if 

I I'd < 1-5, \u v \ < 1.5, and \v^ + v v \ < 1.5 (4.75) 

As depicted in Fig. 14, the domain of stability defined by Eq. (4.75) is a hexagonal 
region in the u^-i ' v plane. Moreover, it will also be shown later that Eq. (4.75) can be 
interpreted as the requirement that the physical domain of dependence of Eq. (4.1) 
be within the numerical domain of dependence. Note that the points on the 
plane that violate Eq. (4.52) form the boundary of a hexagonal region which is entirely 
within the stability domain defined in Eq. (4.75). As was emphasized earlier, the 2D 
a scheme applies even at these points. 

(d) It is shown in [9] that the 2D a scheme has the following property, i.e., for any 
(j,k,n) 6 fi, 

q (j, k,n + 1) — » q(j, k,n) as At — ► 0 (4.76) 

if a x , a,y , w , 6, and h are held constant. The ID d scheme also possesses a similar 
property, i.e., Eq. (2.19) in [2]. The above property usually is not shared by other 
schemes that use a mesh that is staggered in time, e.g., the Lax scheme [52]. 

(e) As in the case of the ID a scheme, the 2D a scheme is also a two-way marching scheme. 
In other words, Eqs. (4.53)-(4.58) can also be used to construct the backward time- 
marching version of the 2D a scheme. More discussions on this subject are given in 
[9]. 

This section is concluded with the following remarks: 

(a) the 2D a scheme is only a special case of the 2D a-p scheme described in [9]. It is a 
solver for the 2D convection-diffusion equation 

du du du ( d 2 u d 2 u\ 

m + dx + “• ay ~ 11 (a* + v j =0 (4 - 

where a x , a y , and p (> 0) are constants. Note that this solver, as in the case of its 
ID counterpart, is unconditionally stable if a x = a y = 0. 

(b) It should be emphasized that, with the aid of Eqs. (4.17)-(4.20), (4.22), and (4.23), 
the 2D d scheme can also be expressed in terms of the marching variables and the 
coefficients tied to the coordinates ( x,y ). In other words, the coordinates ((,7]) are 
introduced solely for the purpose of simplifying the current development. The essence 
of the 2D d scheme , and the schemes to be introduced in the following sections, is 
not dependent on the choice of the coordinates in terms of which these schemes are 

expressed. 
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5. The 2D a-e and a-e-a-/3 Schemes 


The 2D a scheme is non-dissipative and reversible in time. It is well known that a non 
dissipative numerical analogue of Eq. (4.1) generally becomes unstable or highly dispersive 
when it is extended to model the 2D unsteady Euler equations. It is also obvious that a 
scheme that is reversible in time cannot model a physical problem that is irreversible in 
time, e.g., an inviscid flow problem involving shocks. As a result, the 2D a scheme will 
be extended to become the dissipative 2D a-e and a-e-a-f3 scheme before it is extended 
to model the Euler equations. As will be shown, the 2D extensions are carried out in a 
fashion completely parallel to their ID counterparts. 

5.1. The 2D a-e Scheme 

To proceed, note that the CEs for the 2D a-e scheme generally are not those associated 
with the 2D a scheme. Here only a single CE is associated with a mesh point ( j,k,n ) € D. 
This CE, denoted by CE(j,A:,n), is the union of CE r (ji,fc,n), r = 1,2,3. In other words, 

CE(j, k,n) d = [<7Ei(j,fc,n)] U [CE 2 {j,k,n)\ U [CE 3 (j,k,n)\ (5.1) 

Instead of Eq. (4.11), here we assume the less stringent conservation condition 

( £ h* • ds = 0 (5-2) 

JS(CE(j,k,n)) 

Obviously, (i) E 3 can be filled with the new CEs, and (ii) the total flux leaving the boundary 
of any space-time region that is the union of any new CEs will also vanish. 

Moreover, because of Eq. (5.1), Eq. (5.2) must be true if Eq. (4.11) is assumed. As 
a matter of fact, a direct evaluation of Eq. (5.2) reveals that it is equivalent to Eq. (4.65) 
(Eq. (4.68)) if (j',fc,n) € Di (( j,k,n ) £ D 2 ). As a result, Eqs. (4.65) and (4.68) are shared 
by the 2D a scheme and 2D a-e scheme. Recall that Eq. (2.7) is also shared by the ID 
a and a-e schemes. In this section, using a procedure similar to that which was used 
to extend the ID a scheme to become the ID a-e scheme, the two marching steps that 
form the 2D a-e scheme will be constructed by modifying the other equations in the 2D a 
scheme, i.e., Eqs. (4.66), (4.67), (4.69), and (4.70). As a prerequisite, first we shall provide 
a geometric interpretation of the procedure by which the second equation of the ID a 
scheme, i.e., Eq. (2.8), was extended to become the second equation of the ID a-e scheme, 
i.e., Eq. (2.13). 

The key step in extending the ID a scheme to the ID a-e scheme is the construc- 
tion of a central difference approximation of du/dx at the mesh point (j,n). The ap- 
proximation is given as the fraction within the parentheses on the extreme right side of 
Eq. (2.12). Consider a line segment in the x-u space joining the two points (xj_i/2>«j"i/ 2 ) 
and (^j+i/2,«',+i/ 2 )- 11 is obvious that the above central-difference approximation is the 
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value of the slope du/dx of this line segment. In the following modification, instead of con- 
sidering a fine segment in the x-u space joining two points, we begin with the construction 
of a plane in the (-q-u space that intersects three given points. 

To proceed, for any (j, k,n) £ fl 9 , q = 1, 2, let 


u 


t n 


def 



n- 1/2 


? 

(j,k;g,r) 


r = 1,2,3 


(5.3) 


By its definition, u [J k;q r ) is a finite-difference approximation of u at ((j, k;q,r),n). With 
the aid of Eqs. (4.24), (4.27) and (4.28), Eq. (5.3) implies that 


u 


i n 





)] 


n — 1/2 


(5.4) 


For both the case q = 1 (see Fig. 15(a)) and the case q — 2 (see Fig. 15(b)), let P, Q , 
and R be the three points in the (-q-u space with their (i) £- and 77 - coordinates being those 
of the mesh points ((j, k; q,r),n — 1/2), r = 1,2,3, respectively, and (ii) their u-coordinates 
being u'^ k . q r p r — 1*2,3, respectively. It can be shown that the plane in the C-q-u space 
that intersects the above three points is represented by 


where 


and 


= «)",*(< - }*<) + «)", t (<7 - kAr,) + 

(5.5) 

r— 1 

(5.6) 

=' I -1 )’ (“0“*!,,», “ “(Mil.1,) M 

(5.7) 

«)”./= (- 1 )’ («!^) 

(5.8) 


The coordinates of the points O and O c depicted in both Fig. 15(a) and Fig. 15(b) 
are (jA(,kAq,uJ k ) and (jA^,kAq,(u c )J k ), respectively. Here u™ k is evaluated using (i) 
Eq. (4.65) if q = 1 and (ii) Eq. (4.68) if q = 2. Equation (5.5) implies that point O c is on 
the same plane that contains points P, Q , and R. Because generally u 7 - k ^ (w c )” fc , points 
O, P, Q and R generally are not on the same plane. Moreover, for every point on the 
plane represented by Eq. (5.5), 


(iD„ = and (I0 C =K) “’‘ (5 - 9) 

As a result of the above considerations, and the fact that the spatial projection of the mesh 
point (j,k,n) G. flq on the (n — l/2)th time level is the centroid of the triangle formed 
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with the mesh points (( j,k] q,r),n — 1/2), r — 1,2,3, one concludes that ( u )j t k‘>( u c)j,k^ 
and {u c n )^ k are central-difference approximations of u , du/d C, and du/dr /, respectively, at 
the mesh point ( j,k,n ). 

To proceed, for any (j, k,n) 6 fl, let 


= ¥(«<)".» *" d «)?,* =' 


(5.10) 


Then the 2D a-e scheme can be defined as follows: For any (j, fc,n) G Di, we assume 
Eq. (4.65) and n 

(“<)'.» = K + )*.‘ + 2e (“< + “ (5 ' n) 


and 


«)”,* = K + )i,fc + 2e K + - < + )L 


(5.12) 


with the understanding that (u a +)" k and («“+)£* are those defined in Eqs. (4.66) and 
(4.67). On the other hand, for any ( j,k,n ) G Ct 2 , we assume Eqs. (4.68), (5.11) and (5.12) 
with the understanding that (u“ + )£ fc and («$+)£* are those defined in Eqs. (4.69) and 

(4.70). 

With the aid of Eqs. (5.4), (5.7), (5.8), (5.10), (4.66), (4.67), (4.69) and (4.70), it can 
be shown that (i) 


(u'+ -u a +)l k 


, \ n — 1/2 / , 


ra- 1/2 


(5.13) 


and 


« - Ol.t 


, \ n — 1/2 / 

U - 2ut + 4u+ ) - (u-2u+ -2tt+J 


n- 1/2 


(5.14) 


if ( j,k,n ) e fli; and (ii) 

O - Oh = l 


u -f- 2u\ 2 u/( 


\ n — 1/2 / 


u — 4 ui + 2u+ 


\n-l/2 

/ (j,fc;2,2) 


(5.15) 


and 


u c + _ 


/ \ n — l/2 / , \ n 1 / 2 

(. + j.f + a^) -(. + 2.f-4O 0iW) 


(5.16) 


if (j,fc,n) G n 2 . Note that (u ^ )* fe , («f)?, fc , «);,* and K+)” fc are explicitly dependent 
on i/ c and i/, (and therefore exphcitly dependent on At). However, according to Eqs. (5.13)- 
(5 16), (n£ + — u“ + )” fc and (u^ + ~'a° n + ) r j k are ^ ree f rom this depenency. Note that a similar 
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occurrence was encountered in the construction of the ID a-e scheme (see the comment 
given following Eq. (2.14)). 

At this juncture, note that: 

(a) The 2D a-e scheme becomes the 2D a scheme when e = 0. 

(b) For the special case with e = 1/2, Eqs. (5.11) and (5.12) reduce to (u+)? fc = (u c c + )^ k 
and (u+)”- k = (u c +)” k , respectively. 

(c) Using the same reason given in the paragraph preceding Eq. (2.14), one may conclude 
that numerical dissipation in the 2D a-e scheme may be controlled by varying the 
value of e. In fact, it will be shown in Sec. 7 that (i) the 2D a-e scheme is unstable if 
e < 0 or e > 1, and (ii) numerical diffusion indeed increases as e increases, at least in 
the range of 0 < e < 0.7. 

(d) Consider the case ( j,k,n ) £ (lj. Then, with the aid of Eqs. (4.28) and (5.13), 
Eq. (5.11) can be rewritten as: 


(«c) 


71 


6 



) 


71 


+ 


e 

3 


6u 


+ 4 



n—l/2 





2at) 


U r 


) 


n — l/2 -1 


(5.17) 


Let (i) an< ^ ( u v)(j,k/i 2 , 2 ) identified with the values of u, du/dC, 

and du/dr] at the mesh point ((j,k;l,2),n - 1/2), respectively; and (ii) 

( u c)(j,k/i,i ) an< ^ identified with the values of u , du/d( and du/drj at 

the mesh point ((j, k; 1,1), n — 1/2), respectively. Then it can be shown that the ex- 
pression within the brackets on the right side of Eq. (5.17) is 0 ( a <, atj). Furthermore, 
because Eq. (4.26) is applicable only for those points (C,*M) £ SE (j,k,n) only (see 
Figs. 10(b) and 11(b)), the expression enclosed within the first bracket on the right 
side of Eq. (4.26) is 0(A^,Af). From the above considerations, one concludes that 
the error of u*((,rj,t;j,k,n) introduced by adding the extra term involving e on the 
right side of Eq. (5.17) is second order in a(, atj , and At. In other words, addition of 
the term involving e results in lowering the order of accuracy of ( m ^)”*. but not that 
of A similar conclusion is also applicable to Eq. (5.11) for ( j,k,n ) <E fl 2 and to 
Eq. (5.12) for either ( j,k,n ) £ Dj or (j,k,n) £ 0 2 - 

The 2D a-e scheme can also be expressed in the form of Eq. (4.72) if 


/ 


1 - 


Q\ 


(1) def 1 




1 - e 
\ 1-e 
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-(1 - ~ u n )(l + u c ) 

-{l + v c -2e) 

-(1 + - 2e) 
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-(1 -u v )(l + v v )\ 

-(1 +u v - 2e) 

-(1 + u v - 2e) J 


(5.18) 



Q ( 2 ] ^ 


qW d ¥ 


( 1 + V( (l + ^<)(2-^) -(1 + «"c)( 1 + v v)\ 

— (1 — c) —(2 — v<; — 4c) l + i/„-2e 

Vo 0 0 / 

! 1+1^7, — (1 + ^rj)(l + ^c) (1 + l / r;)(2 — I'jj) ^ 

0 0 0 
V _(1 — e) 1 + V( ~ 2e -(2 - v v - 4c) / 


(5.19) 


(5.20) 


/l + U^ + U n (1 + t/( + IS V )(1 - (1 + V{ + ^)(1 "»?)\ 


n (2) def 1 

Vi - 3 


and 


-(i-0 

V -(1-0 


-(1 - U( - 2e ) 
-(1 -u<- 2c) 


-(1 -v v ~ 2c) 

-(1 ~u v - 2c) / 


<#’ = ; 


g', 2) = | 


{X — v^ — (1 — ^)(2 + I'f) (1 ,/ c)(l 

1 - e —(2 + i/£ — 4c) 1 - ^ - 2c 

Vo 0 0 / 

/l — IA? (1 — ^jj)(1 - ^c) — (1 — ^ 

0 0 o 

V 1 — e 1 — i/£ — 2c —(2 + t'jy — 4c) / 


( 5 . 21 ) 


(5.22) 


(5.23) 


Note that, with the above definitions, Eqs. (4.73) and (4.74) are also valid for the 2D a-c 
scheme. 


5.2. The 2D a-e-a-/3 Scheme 

For the same reason that motivates the extension of the ID a-c scheme to become the 
ID a-e-a-P scheme, the 2D a-e scheme will be extended to become the 2D a-e-a-/? scheme. 
As a preliminary for these extensions, first we shall provide a geometric interpretation of 
the procedure by which the ID a-e scheme was extended to become the ID a-e-a-(3 scheme. 

The key step in extending the ID a-e scheme to ID a-e-a-/3 scheme is the construction 
of a nonlinear weighted average of (u^)" and (u^t)” (see Eqs. (2.56)-(2.61)). Let Pj _ — 
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(* C i-l/ 2 ? <U j_i/ 2 )) Pj — i x j^ u J) an d Pj+ — (® j +i /2 5 u j” 1 / 2 ) ke three points in the x-u 
space. Then according to Eqs. (2.12) and (2.56), «t)?, (u c x \)J and ( u '+)?, respectively, 
are equa l to the values of the slope du/dx of the three line segments Pj-Pj, P } Pj + and 
Pj-Pj+i multiplied by the normalization factor ax/4. Equation (2.57) states that ( u c x + )J 
is the simple average of (u c x \)J and (u c x t)]. Thus one can say that the key step in 
extending the ID a-e scheme to become the ID a-c-a-0 scheme is the construction of the 
weighted average of the normalized slopes of Pj-Pj and P } P J+ using the function W D . 
In the construction of the 2D a-e-a-/? scheme, paralleling the evaluation of the values of 
du/dx along the three edges of the triangle AP j _P j P j+ in the x-u space, we shall study 
the gradient vectors Vw associated with the four faces of a tetrahedron in the tj-u space. 
The vertices of the tetrahedron are the points 0, P, Q and R depicted in either Fig. 15(a) 
or Fig. 15(b). The nonlinear weighted average used in the 2D a-e-a-fd will be constructed 
using three of the four gradient vectors referred to above. 

To proceed, consider ( j,k,n ) <E Cl q . Also let planes # 1 , #2, and #3, respectively, be 
the planes containing the following trios of points: (i) points 0, Q, and R ; (ii) points 0, 
R ’ an d P\ and (iii) points 0 , P , and Q. Then; in general, these planes differ from one 
another and from the plane that contains points P , Q and R. In the following derivations, 
first we shall derive the equations representing the former three planes. 

As a preliminary for the developments in this and the following sections, for any real 
numbers s i, s 2 and S 3 , let 


^( 5 l>S2i>S3) — — (2 s 2 + -S3)/a^, 


f\ ^(• s l 1 s 2i s 3) — (2sj +S3 )/a£, 

fl 3 \s 1 ,s 2 ,s 3 ) d = (sj - s 2 )/a(, 


$2,33) = f -(^2 + 2s 3 )/at] 

f q 2 \^l,S 2 ,S 3 ) ^ ( 5 1 - S 3 )/ATJ 

f q 3 \ s 1,82,83) = f ( 2 si + s 2 )/At] 

A'H*,*,;) + s s ), /<■>(„, S 2 , S3 ) * ( 3 * + ”)»» + (3 *-<■>)«» 

//in Iwh 


fL 2) (‘uS„S + 


2 wh 


/< 3 >( S „3 2 ,3 3 ) d =' P-, 4 3 >(3 1>52 ,3 3 ) d =' fy . " 3 *)** +aiP< » 


(5.24) 

(5.25) 

(5.26) 

(5.27) 

(5.28) 

(5.29) 


let 


2 W ’> ' 2 wh 

In the following, consider a mesh point (j, k,n) e fl q (q = 1,2). For any r = 1 , 2 , 3 , 

*r d = (-l)«(«?,Jk - «(,• Jk W , r) ) (5.30) 

( v X ) )?,k = f /< r) («i,* 2 ,® 3 ), K r) )7, fe = f 4 r) («i,x 2 ,x 3 ) (5.31) 

fi'H* («i r) )? tfc = /S r) (* n* 2 ,x 3 ) (5.32) 
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Then it can be shown that, for each r = 1,2,3, plane #r is represented by 


« = (»< r) )”> (C - JA<) + (<’)”,» (•? - *"») 

+ 

if the coordinates (Ci 7 ?) are used; or by 

u = ( u i r) )],k ( x ~ x i,k) + ( ui y r) )],k ( V ~ Vj,k) 
+ «?,* 


(5.33) 


(5.34) 


if the coordinates (x,y) are used. 

Using Eqs. (5.33) and (5.34), one concludes that, at any point on plane #r, r = 1,2,3, 
we have 

/ C% \ / \ 

(5.35) 


and 


9U \ -(Jr))" _ nd (fa\ = ( tt (rhn 
dc) v ~ { < )hk \dr]) c v v h ' J 

= (^>)h (|) =«’)"• 


du 

dx 


(5.36) 


As a result of Eqs. (5.35) and (5.36), at any point on plane # r, r = 1,2,3, (u { x r) )* k and 
i u< y'^)’jk can considered as the covariant components of the vector Vu with respect to 

the Cartesian coordinates (x,y), while (ti ( ( r) )£ fc and ( are the covariant components 
of Vm with respect to the non-Cartesian coordinates (CiV) [55]. Furthermore, according 
to Eq. (5.36), at any point on plane #r, r = 1,2,3, we have 


Vu 


(Or)l k = 


\/ («L r) 


) 2 + (<4r’) 




(5.37) 


Note that, by definition, (0 r )£ fc , r = 1,2,3, are scalars. For readers who are not familiar 
with tensor analysis, it is emphasized that generally (0r)j,fc would not be a scalar and 

therefore the first equality sign in Eq. (5.37) would not be valid if ul * and Uy * in the same 
equation, respectively, are replaced by and u'jK 
To proceed further, let 


/„(r)+xn defAC, (rhn ( (r)+ ^ def ^ / (r)^ 

Then Eqs. (5.7), (5.8), (5.10), (5.24)-(5.26), (5.30) and (5.31) imply that 


(5.38) 


1 


K + )".* = 3 


u™ + +u™ + + u 


( 3 ) + ]” 

< h,k 


(5.39) 


NASA/TM — 1 998-208843 


43 



and 


(5.40) 


« +w 




3 L 


u 


(D+ +u (2)+ 


+ u '% 


i.e., (i) is the simple average of ■u^ r)+ , r = 1,2,3. and (ii) is the simple average of 

tt^ + , r = 1,2,3. Equations (5.39) and (5.40) can be considered as the natural extension 
of Eq. (2.57). Note that, for simplicity, in the above and hereafter we may suppress the 
space-time mesh indices if no confusion could occur. 

Note that, as a result of Eq. (5.38), at any point on plane #r,r = 1, 2, 3, (u^ + )” fc and 

(u^ + )J k are the normalized covariant components of Vu with respect to the coordinates 
(C,v)- On the other hand, as a result of Eqs. (5.9) and (5.10), at any point on the plane 
that contains the triangle A PQR, (“£ + )j,*, and (w^ + )" ife are the normalized covariant 
components of Vm with respect to the same coordinates Recall that planes #1, 

#2, and #3, respectively, are the planes that contain the triangles A OQR, A ORP and 
A OPQ. The last three triangles and A PQR are the four faces of the tetrahedron OPQR. 
Thus Eqs. (5.39) and (5.40) state that Vn associated with one face of this tetrahedron is 
one third of the sum of Vu associated with the other three faces. This conclusion is true 
only because the spatial projection of point O on the plane that contains A PQR is the 
geometric center of A PQR. 

To proceed further, given any a > 0, the nonlinear weighted averages (■u™ + )? fc and 
(tt“ + )" jfc are defined by 


0, 


if $ 1 = e 2 = 0 3 = o 


w+ def 

U, ^ = 


( 0 2 fl 3 ) a t 4 1 )+ + wr + 

(#l# 2 ) a + (#2#3) a + (03#1 )" 


( 2 )+ 


,( 3 )+ 


, otherwise 


(5.41) 


and 


0, 


if $ 1 = e 2 = 0 3 = o 


w-\- def 
= 


(M 3 ) a 4 1)+ + (030i) a < ;+ + (g 102 )°^ 

(^l^ 2 ) a + (Ms)" + { 8381 )“ 


au,( 2 )+ 


)a„< 3 )+ 


, otherwise 


(5.42) 


respectively. To avoid dividing by zero, in practice a small positive number such as 10 -60 
is added to the denominators of the fractions on the right sides of Eqs. (5.41) and (5.42). 
Note that, in the above weighted averages, the weight assigned to a quantity associated 
with plane # r is greater if 8 r is smaller. 

Also note that the above denominators vanish if a > 0, and any two of 8 \ , 62, and 
O3 vanish. Thus, consistency of the above definitions requires proof of the proposition: 
= 82 = 83 = 0, if any two of 8 lf 0 2 , and 83 vanish. 
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Proof: As an example, let 6\ = 02 — 0- Then Eq. (5.37) implies that Ux ^ — u\j ^ — 0, 
r — 1,2. In turn, Eqs. (5.27), (5.28) and (5.32) imply that Xl = x 2 = x 3 - 0. 0 3 - 0 now 
follows from Eqs. (5.29), (5.32) and (5.37). QED. 

As a result of Eq. (5.41), we have 



4 ,,+ . 

if 9i 

= 0, 

02 

> 0, 

and 03 

> 0 

II 

+ 

5* 

«f + , 

if 02 

= 0, 

0i 

>0, 

and 03 

> 0 


4 !)+ , 

if 03 

= 0, 

9i 

>0, 

and 02 

> 0 


Assuming 6 r > 0, r = 1, 2, 3, we have 

(l/*,) a t£ )+ + (l/0 2 ) a u{ 2)+ + (l/^) a « ( C 3)+ 


(5.43) 


(5.44) 


Thus the weight assigned to is proportional to (l/0 r ) a . By using (i) Eqs. (5.39), 

(5.41) and (5.44), and (ii) the fact that = 0, r = 1,2,3, if 9 r = 0, r = 1,2,3, one 

arrives at the conclusion that 


if = 02 = 03 (5.45) 

Obviously Eqs. (5.43)-(5.45) are still valid if each symbol ( is replaced by the symbol q. 

With the above preliminaries, the 2D a-e-a-/? scheme can be defined as follows: For 
any ( j,k,n ) 6 Oi, we assume Eq. (4.65) and 



n 

j 1 * 


= («< + >”.* +2< ( u < + 



(5.46) 


and 


(Oh = Oh + 2e « - Oik + w + - O”. 


(5.47) 


with the understanding that (u£ + )£fc and (u“ + )” k are those defined in Eqs. (4.66) and 
(4.67). On the other hand, for any (j,k,n) £ fi 2 , we assume Eqs. (4.68), (5.46) and (5.47) 
with the understanding that (u“ + )” fc and («“+)”* are those defined in Eqs. (4.69) and 
(4.70). 

At this juncture, note that, on the smooth part of a solution, 9 \ , 0 2 , and 03 are nearly 
equal. Thus the weighted averages and are nearly equal to the simple averages 

and u c v + , respectively (see Eq. (5.45)). As a result, the effect of weighted- averaging 
generally is not discernible on the smooth part of a solution. 

Finally note that, according to Eq. (5.37), evaluation of (0 r ) Q does not involve a 
fractional power if ct is an even integer. Because a fractional power is costly to evaluate, 
use of the a-e-a-/3 scheme is less costly when a is an even integer. 
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6. The Euler Solvers 


We consider a dimensionless form of the 2-D unsteady Euler equations of a perfect 
gas. Let p , u , v, p, and 7 be the mass density, x-velocity component, y-velocity component, 
static pressure, and constant specific heat ratio, respectively. Let 

u i -Pi u 2 = pu, u 3 =pv, u 4 = p /(7 - 1) + p(u 2 + v 2 )/2 (6.1) 


ft = 

ft = (7 - 1)«4 + (3 - 7)(u 2 ) 2 /( 2 «i) - (7 - l){u z ) 2 / {2ui) 

f 3 = u 2 U 3 /U! 

ft = 7 u 2tt4/«l - (l/2)(7 - l)u 2 [(«2) 2 + («3) 2 ] /(ui) 2 

fl = «3 
ff = U 2 Uz/Ui 

ft = (7 - 1)«4 + (3 - 7 )M7(2«,) - (7 - l)(u 2 ) 2 /(2ti,) 

and 

ft = 7u 3 u 4 / Ul - (l/2)( 7 - 1)« 3 [(« 2 ) 2 + («s) 2 ] /M 2 
Then the Euler equations can be expressed as 


du m Offn 9fl = 

dt Ox dy 


m = 1 , 2 , 3, 4 


( 6 . 2 ) 

(6.3) 

(6.4) 

(6.5) 

( 6 . 6 ) 

(6.7) 

( 6 . 8 ) 

(6.9) 


( 6 . 10 ) 


Assuming smoothness of the physical solution, Eq. ( 6 . 10 ) is a result of the more funda- 
mental conservation laws 


(p h m ■ ds = 0, m = l,2,3,4 (6.11) 

Js(V) 

where 

hm = {fmifmi u m) , m = 1,2, 3, 4 (6.12) 

are the space-time mass, i-momentum component, y- momentum component, and energy 
current density vectors, respectively. 

As a preliminary, let 

fm,i d = dfm/due, and f v m /= dfyjdu t , mj= 1,2, 3, 4 (6.13) 

The Jacobian matrices, which are formed by f^ ( and e , m,£ = 1,2, 3, 4 , respectively, 
are given in [ 9 ]. 
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Because /£ and /£, m = 1,2, 3, 4, are homogeneous functions of degree 1 [53] in u u 
U 2 , Ui, and U 4 , we have 


4 •* 

fm = fra, l U (-i & nd ffn ~ 


(6.14) 


Note that Eq. (6.14) is not essential in the development of the CE/SE Euler solvers to be 
described in the following subsections. However, in certain instances, it will be used to 
recast some equations into more convenient forms. 

6.1. The 2D Euler a Scheme 

For any (x,y,t) G SE(j,fc,n), u m (x,y,t), /*(*,»,<), /&(*»»,*)> and h m (x,y,t), re- 
spectively, are approximated by u* m (x,y,t-,j,k,n), f^(x,y,t-,j,k,n), f%?{x,y,t;j,k,n), 
and h* m (x,y,t‘,j,k,n). They will be defined shortly. Let 

U* m (x,y,t',j,k,n) 1= (u m )fk + ( u mx)fk{ X ~ + ( U my)j,k(y ~ ) (6.15) 


+ (timt)j, k(t t n )i 


m — 1 , 2 , 3, 4 


where (u m )? |fc , (»«„)£*, and (u ml )? tfc are constants in SE(j ,*,«)• Obviously, 

they can be considered as the numerical analogues of the values of u m , dum/ox , oum/oy , 
and du m /dt at (*j,fc,2/j,fe,< n ), respectively. 

Let (/*)£*, (/&)£*, (/*,/)£*, and (/*,<)?,* denote the values of /* , /» , /£,<> and 
respectively, when ii m , m = 1,2, 3, 4, respectively, assume the values of (i*m)", fc i 
m = 1,2, 3,4. For any m = 1,2, 3, 4, let 

(6 ' 16) 

< 6 - n > 

*=1 

(/£,)?,/= £(/»,<)", »(“«)".* (6 ' 18) 

(«» = r (fU9) 
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and 


UDlk = £ ( 0 “.*(««) 2 » 


t= 1 


Because (i) 


8fi 


dx 


= Y f x 

/ > J m,i 


1=1 


due 

lh' 


m = 1,2, 3, 4 


( 6 . 21 ) 


(6.22) 


and (ii) the expression on the right side of Eq. (6.16) is the numerical analogue of that 
on the right side of Eq. (6.22) at (x ltk ,y hk ,t n ), can be considered as the nu- 

merical analogue of the value of dfa/dx at (x i)fc ,y iifc ,r). Similarly, (fa y )] k , 

UUlkAfly)l k , and can be considered as the numerical analogues of the values 

°f dfrh/dy, df^/dt, df^/dx, df^/dy, and df-^/dt at (xj, k ,yj, k ,t n ), respectively. As a 
result, we define 


/£•(*, — ( fm)J,k + (fmx)],ki X - X i,k) + (fmy)j,k(y ~ Vj.k) 


+ (/„«)"*(* - n. 


m = 1,2, 3, 4 


(6.23) 


and 


fm( x ,y,i',j,k,n) d = (fl)l k + (fl x )l k (x - x jtk ) + (fl y )l k (y - y hk ) 

+ (fUhit - n, m = 1 , 2 , 3, 4 

Also, as an analogue to Eq. (6.12), we define 

h' m (x,y,t;j,k,n) d = (fZ>‘(x,y,i;j,k,n),fl’(x,y,t;j,k,n), 


(6.24) 


u. 


i(®> y> t ; j, A:, 7i)^ , m — 1,2, 3, 4 


(6.25) 


Note that, by their definitions: (i) (/»)£*, (/»)"*, (/*,<)?,*, and (/» ^ are functions of 
( u m)j ik ,m - 1,2, 3, 4; (ii) (/£*)£* and (/&„)£*. are functions of (um)^ and (u mx )J ik ,m = 
1,2, 3, 4; (iii) (f^ y )^ k and {f* iy )j, k are functions of (tt m )£ fc and (u m y )* k , m = 1,2,3, 4; 
and (iv) (/£ t )£ fc and (f^ t )^ k are functions of (u m )£ fc and ( u mt )J k , m = 1,2, 3, 4. 
Moreover, we assume that, for any ( x,y,t ) € SE(j,A;,n), and any m = 1,2, 3, 4, 


^<(j,t/,f;j,fc,n) + d/**(x,y,*;j,A:,n) + df£ (x, y, t ; j, k, n) = 




dx 


dy 


(6.26) 


Note that Eq. (6.26) is the numerical analogue of Eq. (6.10). With the aid of Eqs. (6.15), 
(6.23), (6.24), (6.16), and (6.20), Eq. (6.26) implies that, for any m = 1,2, 3,4, 

4 

( U mt)j >k = —(fmx)j, k ~ (fmy^k = ~ ^ [/m,< u tx + **<»] . (6.27) 
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Thus {u mt )l k is a function of and (u my )J )fc . From this result and the 

facts stated following Eq. (6.25), one concludes that the only independent discrete variables 
needed to be solved for in the current marching scheme are {u m )J k , {u m3 ,)^ k , and (u my )? k . 

Consider the conservation elements depicted in Figs. 10(a) and 11(a). The Euler 
counterpart to Eq. (4.11) is 


h* m -ds = 0, r = 1,2,3, m = l,2,3,4 


(6.28) 


S(CE T (j,k,n)) 


Next we shall introduce the Euler counterparts of Eqs. (4.22), (4.23), (4.27), and 
(4.28). For any (j,k,n) £ ft, let 


(C,,)h 


def j,-l 


(O?,* 


m, i — 1,2, 3, 4 


(6.29) 


auu . 

/ («mc)j,fc \ / ( u mx)y k \ 

j d = T* j , m — 1,2, 3, 4 (6.: 

\(u mT ))jjt/ \( U my)j y k J 

The normalized counterparts of those parameters defined in Eqs. (6.29) and (6.30) are 


(6.30) 


( f<+ \n \n an J ( P+ \ri fV \n 

(im,f — 2 aC anCl 2 a 


(6.31) 


, 4. , n def ( \n J t + \n def / \n 

(' U m()j,k = “ g-( tt mc)j,fc ’ and ( U m V )j,k ~ U ™v)j,k 


(6.32) 


In the following development, for simplicity, we may strip from every variable in an 
equation its indices j , k , and n if all variables are associated with the same mesh point 
( j,k,n ) £ ft. Let F‘> + and F v+ , respectively, denote the matrices formed by and 

f n+ . m l = 1,2, 3, 4. Let I be the 4 X 4 identity matrix. Then the current counterparts 

J m,£’ ‘ ‘ 7 7 

to Eqs. (4.29)-(4.46) are 



(6.33) 

d = ±(J - F c+ - F 7)+ )(J + F c+ ) 

(6.34) 

d = ±{I - F<+ - F’+KJ + F’» + ) 

(6.35) 

d 4 f / + F<+ 

(6.36) 

E& )=t d = =f(F + F^ + )(2I - F<+) 

(6.37) 
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s (l)± def ± ( T + F C+)(J + F -J+) 


(6.38) 

= f I + F v+ (6.39) 

s£ ):t d = ±(J + F”+)(J + F< + ) (6.40) 

S^ 3 ):t d = q=(J + F’ 7+ )(2/ - F”+) (6.41) 

s (2)± def 7 + F c+ + F ,+ (6.42) 

Eg* d J: f T (J + F<+ + F* + ){I - F<+) (6.43) 

Eg ):t d = =f(F + F<+ + F V+ )(I - F v+ ) (6.44) 

s (2)± def 7 _ F c+ ( 6 .45) 

Sg ):fc d = ±(J - F<+)(2J + F <+ ) (6.46) 

Eg )±d = T(/-F< + )(J-F , ' + ) (6.47) 

s (2)± def J _ F „+ (6.48) 

sS 2 2 }± d = =F(jT - F 7?+ )(J - F<+) (6.49) 

and 

sg ):£ = f ±(I - F"+)(2/ + F^) (6.50) 


Note that Eqs. (4.29)-(4.46) become Eqs. (6.33)— (6.50), respectively, under the following 
substitution rules: 

§1: 1, */f, and u n , be replaced by J, F^ + , and F v+ , respectively. 

§2: cril^ be replaced by q = 1,2 and r, s = 1,2,3, respectively. 

As will be shown, under the above and other rules of substitution to be given later, many 
other equations given in Secs. 4 and 5 can be converted to their Euler counterparts given 
in this section. The latter will be referred to as the Euler images of the former. 

Equation (6.28) is evaluated in Appendix C. Let (j, k,n) £ f l q . Let u, u t , u+, and «+, 
respectively, be the 4x1 column matrices formed by u m , u mt , u + and u m = 1, 2, 3, 4. 
Then, with the aid of Eq. (6.14), for any pair of q and r (q = 1,2 and r = 1,2,3), the 
results with m = 1,2, 3, 4 can be combined into the matrix form 


x[\ )+ e + zTi z c + + s!l )+ «+ 


Eq. (6.51) is the Euler image of Eq. (4.51) under the substitution rules §2 and 
§3: u, u t , u^, and «+ be replaced by u, u t , , and tt + , respectively. 
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As a result of Eqs. (6.33)-(6.50), we have 


v (g)± , v (?)± , y(9>± _ o T 

Zj 11 -h ^21 ' ^31 — 01 > 


9 = 1,2 


(6.52) 


E< 9)± + E< 9)± + E ( 3 9)± = S< 9)± + E?* + E ( 3 9)± =0, 9 = 1,2 (6.53) 

Equations (6.52) and (6.53) are the Euler images of Eqs. (4.47) and (4.48), respectively. 
For either q = 1 or q = 2, by summing over the three equations r = 1,2,3 given in 
Eq. (6.51), and using Eqs. (6.52) and (6.53), one concludes that, for any (j,fc,n) € ft,, 




Et , 1 ) "« + s ( r 9) -tr+ + E ( 4 ) u+ 




, 9 = 1,2 


(6.54) 


As a result, u 7 ?. can be evaluated in terms of the marching variables at the (n - l/2)th 
time level. 

Note that, with the aid of Eqs. (6.33)-(6.50), Eq. (6.54) can be expressed exphcitly as 

* = !((/- + d + + d + **)£?„ 4 °] 

S L (6.54a) 

if ( j,k,n ) E fti; or 

* = 1 [</ + »; 2) + d - * + )r;&> 4 2) + d - ^ ; & *’•] 

(6.546) 

if (j,fc,n) € fl 2 - Here (i) 


3j 1) d = [«-(/ + F<+) « c + - (/ + F ,+ ) 


-+( 1 ) def 


’« + (2/-F <+ )tT+ - (J + F^) 


n— 1 / 2 


-(1) def 


W [«-(/ + F< + ) tt+ + (21 - F v+ ) 


n — 1/2 


with (jf, fc,n) £ fii ; and (ii) 


f > d 4 f [«+(/- F<+) u+ + (/ - F v+ ) u. 


(j,fc;2,l) 


(2/ + F c+ )u+ + (/-F” + ) 


(j>fc;2,2) 
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and 


-(2) def 


U 


n- 1/2 


+ (I - F <+ ) - (2/ + F- + ) ufl (6.60) 

with ( j,k,n ) 6 O 2 . Eqs. (6.54a)-(6.60) are the Euler images of Eqs. (4.65), (4.68) and 
(4.59)-(4.64), respectively, under the substitution rules §1, §3 and 

§4: be replaced by si q \ q = 1,2, and r = 1,2,3, respectively. 

For any (j,k,n) £ the matrices (E^ + )? fc , r = l ? 2, 3, are known functions of 
uj k . Thus they can be evaluted after the latter is evaluated using Eq. (6.54). Assuming 

the existence of the inverse of each of the matrices (E^ + )" fc (see Appendix D.3 for an 
existence theorem), it follows that one can also evaluate Sr (q = 1,2 and r = 1,2,3) 


where 


5 ( 9 ) 4| f 


(s ( r ? ,+ )" 

-1 

X 




1 n- 1/2 




(6.61) 


Note that, in this paper, the inverse of a matrix A is denoted by [A] -1 . 

At this juncture, note that S ^ can be evaluated by a direct application of Eq. (6.61), 
if one does not mind inverting the 4x4 matrices . Alternatively, for each pair 

of q and r, one may use the method of Gaussian elimination to obtain the 4x1 column 


* ( n\ 

matrix Sr } as the solution to the matrix equation 


( E ‘* )+ )m = + E rt~“< + + 

Furthermore, by multiplying Eq. (6.51) from the left with 


-in-l/2 
- (j,k',q,r) 


(6.62) 




i -1 


repeatedly with all possible pairs of q and r, and using Eqs. (6.33)-(6.50) and (6.61), one 
has [9] (i) 

r 1 

f(l) 


u + 


(1) 


and 


(I + F< + )u+ + (I+F r > + )u + 1” = §} 
u- (2I-F< + )u+ + (I + F" + )u+] n =S. 

u + (I + F< + ) u+ - (21 - F v+ ) u +1 ” = 5 3 (1) 
where (j, k,n) £ flj; and (ii) 

(I - i* + ) u+ -(I- F r,+ ) u] 


in = V (2) 
j,k 1 


(6.63) 

(6.64) 

(6.65) 

( 6 . 66 ) 
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and 


u+ (2 I + F <+ )uf 


(I - F v+ ) 


n 



u- (I-F<+)u+ + (2 1 + F^u+r 


_ c(2) 
— *^2 

(6.67) 

C (2) 

— 

(6.68) 


where ( j,k,n ) G f^ 2 - 

Note that, with the aid of Eqs. (6.33), (6.36), (6.39), (6.42), (6.45), (6.48) and (6.61), 
Eq. (6.54) can also be expressed as 


”1* = ^ [(' - F" ~ S , , <,) + (I + ^ + )“, +(I + F'Xk 

if ( j,k,n ) € fii; or 

= 5 [(' + F<+ + F, % s' 2 ’ + < 7 - + (' “ f ’ + )m s ' 


n 5(2) 
3 


(6.69) 


(6.70) 


if ( j,k,n ) G fi 2 . Furthermore, by subtracting Eqs. (6.64) and (6.65), respectively, from 
Eq. (6.63), one obtains 



and 



respectively, where ( j,k,n ) G Dj. 
(6.68), respectively, one obtains 


— (,7“+)™ d — f 
~ )j,k - o 


(X) _ £(») 


) 


(6.71) 


= K + )”„ = l (s, (1) - s 3 (1) ) 


(6.72) 


Next, by subtracting Eq. (6.66) from Eqs. (6.67) and 



(,y a +) n ^ 

v“C lj,fc — 3 ^■- , 2 



(6.73) 


and 




— ii7“+) n 

\ U r/ ) j,k 


1 

3 




(6.74) 


respectively, where (j,k,n) G Cl 2 - 

Note that, under the substitution rules §1, §3, 

§5: u“ + and u“+ be replaced by n“ + and u“ + , respectively. 

§6: si q) be replaced by S r (<?) , q = 1,2, and r = 1,2,3, respectively. 

Eqs. (6.63)-(6.74) are the Euler images of Eqs. (4.53)-(4.58), (4.65), (4.68), (4.66), (4.67), 
(4.69) and (4.70), repectively. 

The 2D Euler a scheme is formed by repeatedly applying the two marching steps 
defined, respectively, by (i) Eqs. (6.54a), (6.71) and (6.72); and (ii) Eqs. (6.54b), (6.73) 
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and (6.74). Note that: (i) because Sr^ can not be evaluated without Uj k being known 
first, one cannot evaluate Uj k using Eqs. (6.69) and (6.70); and (ii) the 2D Euler a scheme 
is a two-way marching scheme in the sense that the conservation conditions Eq. (6.28) can 
also be used to construct its backward time marching version. 

At this juncture, note that the 2D Euler a scheme is greatly simplified by the fact that 
u^ k can be evaluated explicitly in terms of the marching variables at the (n — 1/2 )th time 

levels using Eq. (6.54). As a result, the matrices (E^ + )£*, which are nonlinear functions 
of can be evaluated easily. In other words, nonlinearity of the above matrix functions 
does not pose a difficult problem for the 2D Euler a scheme. 

To explain how Eq. (6.54) arises, note that, because of Eq. (5.1), 

<J) h*^ ■ ds — 0, (j, fc,n) G fl (6.75) 

Js(CE(j y k,n)) 

is the direct result of Eq. (6.28), the basic assumptions of the 2D Euler a scheme. According 
to Eq. (5.1), CE(j,k,n) is the hexagonal cylinder A r B'C* D* E 1 F f ABC DEF depicted in 
Figs. 10(a) and 11(a). Except for the top face A! B* C* D 1 E f F * , the other boundaries of this 
cylinder are the subsets of three solution elements at the (n — l/2)th time level. Thus, for 
any m = 1,2, 3, 4, the flux of leaving CE(j, fc,n) through all the boundaries except the 
top face can be evaluated in terms of the marching variables at the (n — l/2)th time level. 
On the other hand, because the top face is a subset of SE(j, fc,n), the flux leaving there is 
a function of the marching variables associated with the mesh point (j, fc, n). Furthermore, 
because the outward normal to the top face has no spatial component, the total flux of 
h ^ leaving CE(j, fc,n) through the top face is the surface integral of u ^ over the top face. 
Because the center of SE(j,k,n) coincides with the center of the top face , it is easy to 
see that the first-order terms in Eqs. (6.15) do not contribute to the total flux leaving the 
top face. It follows that the total flux leaving the top face is a function of ('U m )j k only. 
As a result of the above considerations, uj k can be determined in terms of the marching 
variables at the (n — l/2)th time level by using Eq. (6.75) only. Equation (6.54) is the 
direct results of Eq. (6.75). 

Because implementation of the 2D Euler a scheme requires, at each mesh point 
(j, k,n) 6 fl, the solution of the three matrix equations (corresponding to r = 1,2,3) 
given in Eq. (6.62), the scheme is referred to as locally implicit [1, p.22]. A simplified and 
completely explicit version of it will be described immediately. 

6.2. The Simplified 2D Euler a Scheme 

Eq. (6.75) is assumed in the 2D Euler a scheme. As a result, Eq. (6.54) is also 
applicable to the new scheme. 

To construct the rest of the simplified scheme, note that, with the aid of Eqs. (6.33)- 
(6.50), a substitution of the approximations 



M ? )+r _1/2 

V rl / (j,k;q,r) 


(6.76) 
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into Eq. (6.61) reveals that 

5 r (9) « J*< 9) , q = 1,2; r = 1,2,3 (6.77) 


where are defined in Eqs. (6.55)-(6.60). 


As 

a result of Eq. (6.77), Eqs 

. (6.71) and 

(6.72) 

can 

be approximated by 




= (< + )”.* 

def 1 / 

_ 3 ‘ 


- 4 ") 

(6.78) 

and 









= K' + )",» 

def 1 / 

“ 3 V 

if* 

- 4”) 

(6.79) 

respectively, where (j,k,n) E Di . 

Similarly, 

Eqs. (f 

i.73) 

and (6.74) can be a] 

pproximated 

by 









~ \ u ( )hh 

def 1 / 

_ 3 ‘ 

V 2 > 

s 2 

\ 

-»T ( ) 

(6.80) 

and 








K + )” t 

_ / 7 7*a' + \ri 

- \ u v 

def 1 / 

- 3 V 

s (2) 

s 3 

-if') 

(6.81) 


respectively, where (j,fc,n) E 1^2- 

Note that Eqs. (6.78)— (6.81) are the Euler images of Eqs. (4.66), (4.67), (4.69) and 
(4.70) under the substitution rules §3, §4 and 

§ 7 : u“ + and «“+ be replaced by u“' + and respectively. 

The first marching step of the simplified 2D Euler a scheme is formed by Eqs. (6.54a), 
(6.78) and (6.79). The second marching step is formed by Eqs. (6.54b), (6.80) and (6.81). 

Moreover, because every (and thus every (u£ + )” fc an< ^ ('“jj + )'j,k {j)k,n) € D) 

can be evaluated without solving a system of equations, the simplified version is compu- 
tationally more efficient than the original scheme. 

6.3. The 2D Euler a-e Scheme 

Eq. (6.75) is assumed in the 2 D Euler a-e scheme. As a result, Eq. (6.54) is also 
applicable to the new scheme. As will be shown shortly, by considering their component 
equations separately, the vector equations that form the rest of the 2 D Euler a-e can be 
developed in a fashion similar to that which was used to develop the 2 D a-e scheme. 

Let (j, k,n) E and consider any m = 1,2, 3, 4. Let («m)[‘j,it;,,r)i K)^i ( u mc)J,ki 
and {u c mi be defined by a set of equations identical to Eqs. (5.3) and ( 5 . 6 )— (5.8) except 
that the symbols u\ u, u t , u c , u ® and u c n in the latter equations are replaced, respectively, 
by the symbols (u' m ), u m , u m t, u c m , u c m( ~ and u c mv in the former equations. Let P m , Q m and 
R m (see Figs. 16(a) and 16(b)) be the three points in the (-rj-u space with (i) their <- and 
77 -coordinates being those of the mesh points ((j, fc; q,r),n — 1 / 2 ), r — 1,2,3, respectively, 
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and (ii) their ^-coordinates being w>* r = 1,2,3, respectively. It can be shown 

that the plane in the C,-rf-u space that intersects the above three points is represented 
by an equation that is identical to Eq. (5.5) except that the symbols u c , u £ and u * in 
Eq. (5.5) are now replaced by u c m , u c m ^ and u c mv , respectively. As a result, for every point 
on the plane referred to above, we have two relations that are identical to those given in 
Eq. (5.9) except that the symbols and u c v in Eq. (5.9) are now replaced by u c m ^ and u c mJ] , 
respectively. Furthermore, let and («£+)£* be defined using an equation that is 

identical to Eq. (5.10) except that the symbols u £, u^ + and u n in the latter equation 

are replaced, respectively, by the symbols tt^, u c m( ., «£+ and u c mri in the former equation. 

Moreover, let u' , u c u', u c ^ + and respectively, denote the 4x1 column 

matrices formed by u' m , u c m , u c m< , u c mj) , and u£+ , m = 1,2, 3,4. Then, with the aid of 
the relation 

Ut = + F^+u+^j (6.82) 

which follows from Eqs. (6.27), (6.29), and (6.30), it becomes evident that we can obtain a 
set of equations that are the Euler images of Eqs. (5.3), (5.4), (5.6)-(5.8), and (5.10)-(5.12) 
under the substitution rules §1, §3, §5 and 

§8: u', u c , u £, u c v , and u c v + be replaced by u', u c , u u u \ ; + and u^ + , respectively. 

Note that the Euler images of Eqs. (5.13)-(5.16) under the substitution rules §3, §5 
and §8 are not valid for the current scheme because (i) (n^ + )”* and (it“ + )” fc are defined 

in terms of Sr q \ q = 1,2, r = 1,2, 3 (see Eqs. (6.71)-(6.74)), while (^“ + )j and (tt$) + )" 
are defined in terms of s[ q) , q == 1,2, r ~ 1,2,3 (see Eqs. (4.66), (4.67), (4.69), and (4.70)); 
and (ii) Sr 9 \ which were defined by Eq. (6.61), axe structually different from s[ q \ which 
were defined by Eqs. (4.59)-(4.64). However, as will be shown shortly, the Euler images of 
Eqs. (5.13)-(5.16) under the substitution rules §3, §7 and §8 do exist. 

For future reference, several key equations associated with the 2D Euler a-e scheme 
will be given explicitly. They are: 


-~*f n def . —r t — — p 

U U,ki g,r) = I U + — 


and 


At 

T 


n-l/2 


u — 2 ” 1/2 

L V C ^ / \ (j,k\q,r) 


(ii e+ \ n def ( l) 1 * / ->>n _ yf n \ 

)j,k g \ U (j,k;q, 2) u (j,k;g,l) J 

«>** = (*&»»> - *«*,..>) 
(ff+)" t = («;+)?,, + 2 £ («-+ - 

K + )“* = K + )",* + 2^(“T - 


(6.83) 

(6.84) 

(6.85) 

( 6 . 86 ) 

(6.87) 
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where ( j,k,n ) £ il q , q = 1,2. The 2D Euler a-e scheme is formed by Eqs. (6.54), (6.86) 
and (6.87) for any (j,k,n) £ il q . 

6.4. The Simplified 2D Euler a-e Scheme 

The defining equations of the simplified 2D Euler a-e scheme are identical to those of 
the 2D Euler a-e scheme except that Eqs. (6.86) and (6.87) should be replaced by 


and 


(“< + )” 4 = K' + )m + Hi ? ~ 
K)J» = « '*)lk + 2«K + - 


( 6 . 88 ) 


(6.89) 


respectively. 

Moreover, with the aid of Eqs. (6.78)— (6.81) and (6.83)— (6.85), it can be shown that 

0) 


(Sf - = - 


/ V n— 1/2 / , \ n — !/ 2 

(u + 4u+ - 2u ■+ ) - (u - 2 u+ - 2 u + ) 


and 

f T r c + _ u a ' + )”. = - 

\TJ U 7J Jj,k g 

if ( j,k,n ) £ Hi] and (ii) 


\ n — 1/2 / 

2u + + 4u + ) - (u - 2u + - 2u + ) 


n — 1 /2 




(“T - - s [ ( ff+M < + + 2 <)^, - (•-■ - 4 “< + + 2 < 


i-l/2 


and 




, \ n — 1/2 / . xn - 1/2 

( u + 2u + 2u+ ) — ( u + 2«^ - 4tt J 

V *> n /(j,fc;2,l) V ^ /(j,fe;2,3) 


(6.90) 


(6.91) 


(6.92) 


(6.93) 


if ( j,k,n ) £ D 2 . 

Note that, under the substitution rules §3, §7 and §8, Eqs. (6.90)-(6.93) are the Euler 
images of Eqs. (5.13)— (5.16), respectively. Also note that ( u (; + )j,k> ( u <; ( u v )j,fc an< ^ 

(tt“' + )"t are explicitly dependent on F^ + and F v+ (and, as a result of Eq. (6.31), also 
explicitly dependent on At). However, according to Eqs. (6.90)-(6.93), (u£ — u “ )£*. and 

( u c + - ■u“' + )£ fc are free from this depenency. 

6.5. The 2D Euler a-e-a-(3 Scheme 
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In this subsection, the techniques used in constructing the ID Euler a-e-a-(3 scheme 
and the 2 D a-e-a-fl scheme will be combined and used to construct the 2 D Euler a-e-a-f3 
scheme. 

To proceed, for any (j, k,n) 6 Cl q , any m = 1,2, 3, 4, and any r = 1,2,3, let 

Xm,r = (“I) 9 ~ (“m)(i,it;,,r)] (6.94) 

K()j,t ~ f\ ^( x m,l , *m,2, * 771 , 3 ), ( u ml))J,k ~ f^ r \ x m,l , *m,2, * 771 , 3 ) (6.95) 

( u mx)j,k = fx ^ (*”1,1 »*m,2 5 * 771 , 3 )) ( u iny)J,k = fy r \ x m,l , * 771 , 2 ) *m,3 ) (6.96) 

where f^ r \ f\ r \ fi r \ and are the functions defined in Eqs. (5.24)-(5.29). Note that 
Eqs. (6.94)-(6.96) are the Euler counterparts of Eqs. (5.30)-(5.32), respectively. 

To proceed further, for either ( j,k,n ) 6 Di or (j, k,n) £ D 2 , consider any fixed value 
of m = 1,2, 3 , 4. Let P m , Q m and R m be the three points defined in Sec. 6.3. Let 
Om (see Figs. 16(a) and 16(b)) denote the point in the £-77 -u space with the coordinates 
O^C) ^ 77 , (u m )^k )■ Let planes $ 1 , $ 2 , and ^3, respectively, be the planes containing the 
following trios of points: (i) points O m , Q m , and R m ; (ii) points 0 m , R m , and P m ; and 
(iii) points O m , P m , and Q m . Then it can be shown that, for each r = 1,2,3, plane is 
represented by an equation that is identical to Eq. (5.33) except that the symbols iA r) , u[ r \ 

and u on the right side of Eq. (5.33) are now replaced by and (u m ), respectively. 

Alternatively, the plane can be represented by another equation that is identical to 
Eq. (5.34) except that the symbols ui r \ u[ r \ and u on the right side of Eq. (5.34) are now 
replaced by u mx , u my , and (w m ), respectively. As a result, for every point on the plane 
#r, we have a set of relations that are identical to those given in Eqs. (5.35) and (5.36) 
except that the symbols u^ r \ ui[\ u£\ and u y r ^ in the latter equations are now replaced 

by n mr) 5 and Umy ? respectively. It follows that, at any point on plane ^ 7 *, we 

have 


|V«| = {6 mr )lk = [\/(»K) ! + (u ( 4) 




(6.97) 


Furthermore, let 


' m C ' j<k c \ U m,Ti )j,k ^ ( u mt/ )j,l 


(6.98) 


Then Eqs. (6.84), (6.85), (5.24)-(5.26) and (6.94)-(6.96) imply that 


_ 1 r (1)+ , (2)+ , (3)+l n 

' g [^TnC ' U m£ ' U m(, \ j ^ 


(6.99) 
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and 


u 


c + \n 

mr))jyk 


U 


( 1 ) + 
rnr) 


+ U 


( 2 ) + 
mrf 


+ u 


] n 

i.* 


( 6 . 100 ) 


i.e., (i) u ^ is the simple average of u^ + , r = 1,2,3; and (ii) is the simple average 
of r = 1,2,3. Equations (6.97)-(6.100) are the Euler counterparts of Eqs. (5.37)- 

(5.40), respectively. 

With the above preliminaries, it becomes obvious that and , respectively the 
present counterparts of the weighted averages ■u™ 4 ' and u™ + defined in Eqs. (5.41) and 
(5.42), should be defined by 


iu+ def 
U mC = 


0, 

(0m2 0m3) Q u Lc + )° + (&ml ^m2 )° 

(^ml^m2) a + (0m 2^3)° + 


and 


if #ml = #m2 = #m3 = 0 


otherwise 


(6.101) 


0, 


if 0m\ = ft 


m2 


^ m 


17171 1 (^m2 ^m.3 ) a (^m3^ml ) Q ^rn.^ ~f(^ml^m2) u 


( 2 )+ 


,( 3 )+ 

mr} 


(^ml @m2 ) a + (^2^3)“ + (0m30ml) C 


otherwise 


0 m3 = 0 


( 6 . 102 ) 


respectively. Note that, to avoid dividing by zero, in practice a small positive number such 
as 10~ 60 is added to the denominators in Eqs. (6.101) and (6.102). 

Let u™ + (u™ + ) be the column matrix formed by («*+), m = 1,2, 3, 4. Then, for 
any (j, ifc,n) G O, the 2D Euler a-e-a-(3 scheme is defined by Eq. (6.54) and 



,-C + 




- a i + '>h 


(6.103) 


and 

K + )” k = K + )?,, + 2e« + - )l, +3W + ~ ^ + )”,r (6.104) 

where c and /? are adjustable parameters. 

6.6. The Simplified 2D Euler a-e-a-fl Scheme 

For any (j,k,n) E ft, the simplified 2D Euler a-e-a-& scheme is formed by Eq. (6.54) 

=(^'+)J, i +2e(r < + -< + )4 k+ l3(u < »+-uJ + )” l (6.105) 

V s / J,k 


and 


and 


(*,% = (<+)?,* + 2*K + - <+)?,* + /?« + - ^ + )h (6-106) 
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where e and 3 are adjustable parameters. 

6.7. The 2D CE/SE Shock-Capturing Scheme 

Let e = 1/2 and (3 = 1. Then the 2D Euler a-e-a-fi scheme and the simplified 2D Euler 
a-e-a-(3 scheme reduce to the same scheme. For any (j, k,n) G fl, the reduced scheme is 
formed by Eq. (6.54) and 

(£+)” = (u” +)?,» (6.107) 

and 

(*,% = (tr)h ( 6 - 108 ) 

The above scheme is one of the simplest among the 2D Euler solvers known to the authors. 
The value of a is the only adustable parameter allowed in this scheme. Because this scheme 
is the 2D counterpart of the ID CE/SE shock-capturing scheme and shares with the latter 
all the distinctive features described in Sec. 2.8, it will be referred to as the 2D CE/SE 
shock-capturing scheme. 
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7. Stability 


In this section, stability of the 2D a and a-e schemes will be studied using the von 
Neumann analysis. Note that Eqs. (4.73) and (4.74) are valid for these two schemes if 
the matrices Q ( r q) (q = 1,2 and r = 1,2,3) are defined using Eqs. (5.18)-(5.23) with the 
understanding that e = 0 should be assumed for the 2D a scheme. 

To proceed, let 

M™(0 Ct 0 v ) = + Q(») e (i/»(-w<+®,) + Q W e W *)('<-***) (7.1) 


and 

M ( ' 2 \0^,O v ) d = Q[ 2) e - (i/3)i6<+e ' ) + Q {2) e- (l/3)i ~ 20<+eTl) + q\ 2) e~ ii/3){0 <~ 20,,) (7.2) 

Furthermore, for all (j,k,n) €E Cl, let 

q(j,k,n) = q*{n,6 o 6r,)e iU6(+k0 '\ (* =* V^l, -v < 0 <t 0 n <ir) (7.3) 

where q*(n,0^,0 v ) is a 3 X 1 column matrix (see Sec. 4 in [1]). Substituting Eq. (7.3) into 
Eqs. (4.73) and (4.74), one concludes that: (i) 


q*(n + m,6(,0 v ) 


1 m 


(7.4) 


where n — il/2,db3/2,±5/2,..., and m — 0,1,2,...; and (ii) 

q*(n + m,6 ( ,$ v )= [M^(0 ( ,0 v )M^(0O0v )\ " *,) ( 7 - 5 ) 

where n = 0,±1,±2,..., and m = 0,1,2,.... Equation (7.4) implies that the am- 
plification matrix among the half-integer time levels is M (1) (0 C , 0 r) )M (2) (0^,0 ri )\ while 
Eq. (7.5) implies that the amplification matrix among the whole-integer time levels is 

M( 2 \e o 0 v )M^\e o e v ). 

Let A and B be two arbitrary n x n matrices. Then AB and BA have the same 
eigenvalues, counting multiplicity [54, p.53]. Thus the 3x3 amplification matrix among 
the half-integer time levels and that among the whole-integer time levels have the same 
eigenvalues. These eigenvalues may be referred to as the amplification factors. The ampli- 
fication factors are functions of phase angles 0^ and 0^. In addition, they are functions of a 
set of coefficients that are dependent on the physical properties and the mesh parameters. 
These coefficients are (i) u c and u n for the 2D a scheme; and (ii) u c , u v , and e for the 2D 
a-e scheme. Let A a , A 2 , and A 3 denote the amplification factors. In the current stability 
analysis, a scheme is said to be stable in a domain of the above coefficients if, for all values 
of the coefficients belonging to this domain, and all 0( and 6 , ? with 7T <L 0^,0 , ? 7r , 

|Ai| < 1, |A 2 | < 1, and |A 3 | < 1 (7.6) 
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Consider the 2D a scheme. By using its two-way marching nature and the fact that 
its stencil is invariant under space-time inversion, it is shown in [9] that, for any given 
v v , and 0 V , 

|A a A 2 A 3 | - 1 (7.7) 

It follows from Eqs. (7.6) and (7.7) that the 2D a scheme must be neutrally stable, i.e., 

|Ai| = |A 2 | = | A 3 1 = 1, -7r < d{ y d v < 7 r (7.8) 

if it is stable. In other words, the 2D a scheme is non-dissipative if it is stable. Moreover, 
a systematic numerical evaluation of Aj, A 2 , and A3, for different values of and 

0 V , has confirmed that the 2D a scheme is indeed neutrally stable in the stability domain 
defined by Eq. (4.75). In the following, we shall discuss the meaning of this stability 
domain. 

Let (jf, k,n) 6 D. According to Eqs. (4.73) and (4.74), the marching variables at 
the mesh point (j, fc,n) are completely determined by those of seven mesh points at the 
(n — l)th time level (i.e., the mesh point (j, fc,n — 1), and points A, 5, (7, D, E and F 
shown in Figs. 17(a) and 17(b)). As a result, in this paper, the interior and boundary of 
the hexagon ABC DEF shall be referred to as the numerical domain of dependence of the 
mesh point (j,h,n) at the ( n — l)th time level . Note that the dashed lines depicted in 
Figs. 17(a) and 17(b) are the spatial projections of boundaries of CEs. 

The 2D a scheme is designed to solve Eq. (4.1). For Eq. (4.1), the value of u is 
a constant along a characteristic line. The characteristic line passing through the mesh 
point (j,fc,n) will intersect a point on the plane t — t n-1 . The point of intersection, 
referred to as the backward characteristic projection of the mesh point (j, h,n) at the 
(n — l)th time level, is the “domain” of dependence at the (n — l)th time level for the 
value of u at the mesh point (j,fc,n). It is shown in Appendix D.l that the backward 
characteristic projection is in the interior of the numerical domain of dependence if and 
only if Eq. (4.75) is satisfied. 

At this juncture, note that the concept of characteristics was never used in the design 
of the 2D a scheme. Nevertheless, its stability condition is completely consistent with the 
general stability requirement of an explicit solver of a hyperbolic equation, i.e., the analytic 
domain of dependence be a subset of the numerical domain of dependence. 

Next we consider the stability of the 2D a-e scheme. Recall that the 1-D a-e scheme 
is not stable for any Courant number v if e < 0, or e > 1 [2]. Similarly, the results of 
numerical experiments indicate that the 2D a-e scheme is not stable in any domain on the 
pl ane if c < 0 or e > 1 . For any t with 0 < e < 1 , the 2D a-e scheme has a stability 
domain on the plane. The stability domains for several values of e were obtained 

numerically. As shown in Figs. 18(a)-(c), these domains (shaded areas) vary only slightly 
in shape and size from that depicted in Fig. 14. They become smaller in size as e increases. 

Given any pair of and u v , Aj, A 2 and A 3 are functions of 6^ and 9 V . Let (i) 

|A 3 1 < |A 2 | < |Aj | <0, -7 r < 0 O 9 V < 7 r (7.9) 
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and (ii) Ai = 1 when 6^ = B v = 0. Then Aj can be referred to as the principal amplification 
factor; while A 2 and A 3 are referred to as the spurious amplification factors [1]. In gen- 
eral, the principal amplification factor is the deciding factor in determining the accuracy 
of computations [1]. Specifically, numerical solutions may suffer annihilations of sharply 
different degrees at different locations and different frequencies if numerical diffusion asso- 
ciated with Aj varies greatly with respect to 6 ^ , t'o and u v [7, p.20]. Moreover, note 
that (1 — |A r |) is a measure of the numerical diffusion associated with A r , r = 1,2,3. For a 
given e, let D(e) denote the stability domain of the 2D a-e scheme on the u^-v v plane. Let 


Xr(«) d = max (1 - |A r |), r = 1,2,3; USe 

-7T<0 C ,0 V <ir, 


0 < € < 1 


(7.10) 


Then, for a given e and each r, (1 - |A r |) is bounded uniformly from above by Xr(e)- The 
numerically estimated values of Xr( £ ) are plotted in Fig. 19. From this figure, one con- 
cludes that the numerical diffusion, particularly that associated with A 1? can be bounded 
uniformly from above by an arbitrary small number by choosing an e small enough. Note 
that this property is also shared by the 1-D a-e scheme (see Eq. (3.19) in [2]). Moreover, 
the results shown in Fig. 19 indicate that X 2 ( e ) and X 3 ( c ) are much larger than Xi( e ) * n 
the range of 0 < e < 0.5. Thus, in this range, the spurious part of a numerical solution is 
annihilated much faster than the principal part. Also it is seen that the numerical diffusion 
associated with the principal solution, measured by Xi( € )> increases with e in the range of 
0 < e < 0.7. 

Because of the appearance of nonlinear weighted-average terms in its defining equa- 
tions, stability of the 2D a-e-a-(3 scheme is difficult to study analytically. However, results 
from numerical experiments indicate that the stability domain of this scheme is slightly 
larger than that of the 2D a-e scheme when a > 0 and (3 > 0. 

Before we proceed further, several concepts related to stability need to be clarified. 
First note that, to define a numerical problem, one must specify (i) the main scheme (such 
as any solver described in Secs. 4-6) used in the updating of the marching variables at 
the interior mesh points, and (ii) the auxiliary discrete initial/boundary conditions. Thus, 
generally stability is not a concept involving only the main scheme. 

Next note that use of the von Neumann stability analysis can be rigorously justified 
only if the numerical problem under consideration satisfies a set of strict conditions [1]. 
They include (i) the mesh used should be uniform in both spatial and temporal directions, 
(ii) the main scheme used should be linear in the discrete variables, and (iii) the boundary 
conditions used should be periodic in nature. Because (i) the stability conditions generated 
using the von Neumann analysis are expressed in terms of the coefficients of the discrete 
variables and the mesh parameters only, and (ii) the above coefficients and mesh param- 
eters are constant and independent of the initial/boundary conditions, the stability of a 
numerical problem that satisfies the above strict conditions (i)-(iii) is completely indepen- 
dent of the initial/boundary conditions. For this special numerical problem, stability can 
be considered as a concept involving only the main scheme. 
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For a uniform-mesh linear problem with non-periodic boundary conditions, the stabil- 
ity conditions generated from the von Neumann analysis generally are necessary but not 
sufficient conditions for stability. For such a problem, the initial/boundary conditions may 
have an impact on stability and numerical diffusion. Note that the results given earlier in 
this section are obtained without considering this impact. 

Generally, stability of a nonlinear problem is highly dependent on the initial/boundary 
conditions, and therefore highly problem-dependent. As a result, a discussion of the sta- 
bility of nonlinear solvers without specifying the exact initial/boundary conditions, such 
as that to be given immediately, is inherently imprecise in nature. 

To proceed, for each mesh point (j,fc,n) 6 ft, a local Euler CFL number u e > 0 is 
introduced in Appendix D.2 (see Eqs. (D.32)-(D.35)). This number has the following prop- 
erty: For the flow variables at the mesh point (j, fc,n), its analytical domain of dependence 
at the (n — l)th time level lies within the corresponding numerical domain of dependence 
if and only if u € < 1. According to the results of numerical experiments, both the 2D 
Euler a scheme and the simplified 2D Euler a scheme are generally unstable. However 
the former is only marginally unstable when v e ,max < 1 where Ve^ax the maximum 
value of v € ever reached in a numerical experiment. As a matter of fact, in simulating 
smooth flows, its round-off error often never reaches an appreciable level before the end of 
the simulation run. As for the other solvers described in Sec. 6, stability generally can be 
realized if < 1 and 0*05 < c < 1. However, for a nonsmooth flow problem, stricter 

stability conditions such as v ey max < 2/3, 0.1 < c < 1 and a > 1 may apply. 
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8. Conclusions and Discussions 


The space-time CE/SE method was conceived from a global CFD perspective and 
designed to avoid the limitations of the traditional methods. It was built from ground zero 
with a foundation which is solid in physics and yet mathematically simple enough that one 
can build from it a coherent, robust, efficient and accurate CFD numerical framework. A 
clear and thorough discussion of these basic motivating ideas was given in Sec. 1. 

The ID CE/SE schemes [2] were reformulated in Sec. 2 such that the reader can see 
more clearly the structural similarity between the solvers of the ID convection equation 
Eq. (1.1) and those of the ID Euler equations. In addition, this reformulation also paves 
the way for the construction of the 2D CE/SE schemes and makes it easier for the reader 
to appreciate the consistency between the construction of the ID CE/SE schemes and that 
of the 2D schemes. 

It was shown in Sec. 3 that the basic building blocks of the spatial meshes used in 
the 2D CE/SE schemes are triangles. As a result, these schemes are compatible with the 
simplest unstructured meshes, and therefore are applicable to 2D problems with complex 
geometries. Furthermore, because they are constructed without using the dimensional- 
splitting approach, these schemes are genuinely multidimensional. 

The 2D a scheme, a nondissipative solver for the 2D convection equation Eq. (4.1), 
was constructed in Sec. 4. It is a natural extension of the ID a scheme and shares with 
the latter several nontraditional features which are listed following Eq. (4.74). 

Because a nonlinear extension of a nondissipative linear solver generally is unstable 
or highly dispersive, the 2D a scheme was modified in Sec. 5 to become the dissipative 2D 
a-e and a-e-a-0 schemes before it was extended to model the 2D Euler equations. It was 
clearly explained in Sec. 5 that these 2D dissipative schemes are the natural extensions of 
the ID a-e and a-e-a-/3 schemes, respectively. Moreover, as in the case of the latter schemes, 
numerical dissipation introduced in the former schemes is controlled by the parameters e, 
a and /3. 

A family of solvers for the 2D Euler equations were constructed in Sec. 6. Not only 
are these solvers the natural extensions of the ID CE/SE Euler solvers, but their algebraic 
structures are strikingly similar to those of the 2D a, a-e and a-e-a-(3 schemes. 

Next, stability of the 2D solvers described in Sec. 4-6 was discussed in Sec. 7. It was 
shown that the 2D a scheme is nondissipative in the stability domain defined by Eq. (4.75). 
It was also shown that the necessary stability conditions for the 2D solvers include: (i) 
the local CFL number < 1 at every mesh point, and (ii) l>c>0, a>0 and (3 > 0 if 
applicable. Note that these conditions are also necessary stability conditions for the ID 
CE/SE solvers. 

A summary of the key results of the present paper has been given. It is seen that 
each of the present 2D schemes is constructed in a very simple and consistent manner as 
the natural extension of its ID counterpart. This is made possible because of the present 
development’s strict adherence to its two basic beliefs which were stated in Sec. 1. 
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To evaluate the accuracy and robustness of the CE/SE schemes, the two simplest 
schemes among them, i.e., the ID and 2D CE/SE shock-capturing schemes, will be used 
in Part II [3] to simulate flows involving phenomena such as shock waves, contact discon- 
tinuities, expansion waves and their interactions. The numerical results, when compared 
with experimental data, exact solutions or numerical solutions by other methods, indicate 
that these schemes can consistently resolve shock and contact discontinuities with high 
accuracy. Note that other CE/SE schemes described in this paper have also been shown to 
be accurate solvers for other applications [11,13-17,20,24,26-28]. Furthermore, using the 
present method, Yu et a 1. have successfully constructed several accurate solvers for ID 
and 2D problems with stiff source terms [21,22,32]. 

Note that the ID CE/SE schemes have been extended to become accurate 2D and 3D 
solvers by others without using the current approach. After constructing their ID CE/SE 
solver for the Saint Venant equations, Molls et aJ. [29] construct the 2D version using 
the Strang’s splitting technique [56]. Furthermore, several 2D and 3D non-splitting Euler 
solvers have also been constructed by Zhang et aJ. [57-61] without using triangular or 
tetrahedral meshes. 

The triangles depicted in Fig. 5 are obtained by sectioning each parallelogram depicted 
in the same figure into two triangles. The 2D CE/SE solvers can also be constructed using 
the triangles that are obtained by sectioning each parallelogram into four triangles. These 
solvers along with other CE/SE solvers with non uniform spatial meshes [4] will be described 
in future papers. 

This paper is concluded with a discussion of several other extensions. 

8.1. A sketch of a 3D Euler solver 

The CE/SE method can be extended to three spatial dimensions using the same 
procedure that was used in extending the method from one spatial dimension to two spatial 
dimensions. In the 3D case, at each mesh point, the mesh values of any physical variable 
and its three spatial gradient components are considered as independent variables. Because 
there are four independent discrete variables per physical variable (or per conservation law 
to be solved), construction of the 3D a scheme and the 3D Euler a scheme demands that 
four CEs be defined at each mesh point. As will be shown immediately, this requirement 
can be met by using tetrahedrons as the basic building blocks of the 3D spatial mesh. 

To pave the way, consider the 2D case and Figs. 5 and 6(a). The quadrilaterals GFAB, 
GBCD and GDEF are the spatial projections of the CEs associated with the point G ' . 
The CEs in the 3D case can be constructed in a similar fashion. Consider the tetrahedrons 
ABCD and ABCP depicted in Fig. 20. Points G and H are the centroids of ABC D 
and ABCP , respectively. The two tetrahedrons share the face ABC. The polyhedron 
GABCH is then defined as the spatial projection of a CE associated with a space-time 
mesh point G ' . The CE is thus a right cylinder in space- time, with GABCH as its spatial 
base. The point G is the spatial projection of point G' . 

In a similar fashion, three additional CEs associated with the mesh point G' can be 
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constructed by considering in turn the three tetrahedrons that share with ABC D one of 
its other three faces. 

Note that a structured 3D spatial mesh can be constructed from the tetrahedrons that 
are obtained by sectioning the parallelepipeds occupying a spatial region. The details will 
be given in a separate paper. 

8.2. Concept of Dual Space-Time Meshes and Its Applications 

The mesh depicted in Fig. 4(a) is staggered in time, i.e., the mesh points that have 
the same spatial locations appear only at alternating time levels. In Fig. 21(a), the mesh 
depicted in Fig. 4(a) (referred to as the mesh 1) is superimposed on another staggered 
mesh (referred to as the mesh 2), with the mesh points of the latter being marked by solid 
triangular symbols. The combination of the meshes 1 and 2 shall be referred to as the dual 
mesh. As shown in Fig. 21(b), a CE of a mesh point marked by a triangle may coincide 
with a CE of another mesh point marked by a dot. 

Obviously the ID a scheme can also be constructed using mesh 2. As a matter 
of fact, one can even combine two independent ID a schemes, one constructed on the 
mesh 1, and the other on the mesh 2, into a “single” scheme referred to as the ID dual a 
scheme. Similarly one can also construct the dual ID a-e and a-e-a-/3 schemes. Each of the 
new schemes has two completely decoupled solutions. Without considering this decoupled 
nature in the von Neumann analysis, it can be shown that the resulting amplification 
factors of the dual ID a scheme are identical to those of the Leapfrog scheme as given in 
[52, p.100]. Note that the deficiency of the standard practice that the amplification factors 
of the Leapfrog scheme are obtained without taking into account the decoupled nature of 
its solutions was addressed in Sec. 1. 

Let (j,n) be a mesh point of mesh 1 (mesh 2). Then ( j ± 1/2, n) are mesh points of 
mesh 2 (mesh 1). Recall that uV” 1/2 (see Eq. (2.10)) are defined in terms of the marching 
variables at ( j ± 1/2, n - 1/2), which are on the same mesh with (j,n). Thus the two 
solutions on meshes 1 and 2 of either the dual ID a-e scheme or the dual ID a-e-a-(3 
scheme are decoupled. However, by replacing u 'j±i /2 u7 j±i /2 (which are evaluated 

using Eq. (2.8) with the understanding that j be replaced by j ±1/2) in their construction, 
each of the above two schemes will turn into a new scheme in which the solutions on meshes 
1 and 2 become coupled. The coupling results from the fact that u” and are no ^ 

associated with the same mesh. Note that the solutions of the new schemes generally are 
indistinguishable from (or only slightly more diffusive than) those of the original schemes. 

In [12,25], two implicit schemes for solving the convection-diffusion equation Eq. (1.2) 
were constructed using a dual space-time mesh. In the case that fi = 0, both the above 
implicit schemes reduce to the explicit non-dissipative dual a scheme. As a result, the 
amplification factors of these schemes reduce to those of the Leapfrog scheme if (i = 0. 
Furthermore, these two implicit schemes have the property that their numerical dissipa- 
tion approaches zero as the physical dissipation approaches zero. The significance of this 
property was discussed in Sec. 1. 
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In case that /i > 0, both the above implicit schemes are truly implicit. This implicit 
nature is consistent with the fact that, for /i > 0, the value of a solution to Eq. (1.2) at 
any point ( x,i ) depends on the initial data and all the boundary data up to the time t. 
In other words, generally an implicit scheme should be used to solve an initial/boundary- 
value problem, such as one involving Eq. (1.2) with fi > 0. This requirement becomes 
more important as the diffusion term in Eq. (1.2) becomes more dominant. 

In addition, for both the above implicit schemes, the solution at the mesh points 
marked by dots, through the diffusion term in Eq. (1.2), is coupled with that at the mesh 
points marked by triangles if /x > 0. Also it was shown in [12,25] that, in the pure diffusion 
case (i.e., when a = 0), the principal amplification factors of both the above implicit 
schemes reduce to the amplification factor of the Crank- Nicolson scheme [52]. Note that 
the latter has only one amplification factor. 

The concept of dual space-time meshes also is applicable to the 2D and 3D cases. 
As an example, consider a 2D mesh (the mesh 1) with the mesh points marked by circles 
in Fig. 6(a)-(c). For this case, the mesh points of the mesh 2 are points G, C", E\ G " , 
I" and K". In general, if ( j,k,n ) represents a mesh point of the mesh 1, then ( j,k,n' ) 
represents a mesh point of the mesh 2 if and only if (n — n') is a half-integer. Note that a 
more complete discussion of the concept of dual meshes will be given in Part II [3]. 

Note that not only can the concept of dual meshes be used to construct implicit 
schemes, but it can also be used to implement reflecting boundary conditions (see the 
following paper [3]). In addition, this concept is indispensable in the development of a 2D 
triangular unstructured-mesh CE/SE scheme [31]. 

8.3. A discussion on locally adjustable numerical dissipation 

Consider the ID a-e-a-0 scheme, i.e., the scheme defined by Eqs. (2.7) and (2.60). 
With e, a and (3 being held constant, generally numerical dissipation associated with 
this scheme increases as the Courant number v decreases. To compensate for this effect, 
Eq. (2.60) may be replaced by 

W)? = M + )” + 2<MK + - <+)” +/SMW+ - u'+)] (8.1) 

where e(u) and (3(u) are monotonically decreasing functions of v with e(0) = /3(0) = 0. The 
optimal forms of these functions generally are problem-dependent. The scheme defined by 
Eqs. (2.7) and (8.1) has the property that 

« + )? as At * 0 (8.2) 

With the aid of Eq. (8.2), it is easy to see that the new scheme shares with the a scheme 
the same property Eq. (2.19) in [2], i.e., 

u n+1 -» u ? and (ut)] +1 -+ («+)? as At -> 0 (8.3) 

In the new scheme introduced above, numerical dissipation is controlled by the pa- 
rameters f3(v) and a with the first two being the functions of the convection speed 
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a, the mesh interval Ax and the time-step size At. In similar extensions involving solvers 
of more complicated nonlinear equations, the values of these parameters may vary with 
space and time, and their local values generally will be functions of local values of dynamic 
variables, mesh intervals and time-step size. 
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Appendix A. A CE/SE Solver for the Sod’s Shock Tube Problem 
with Non- Reflecting Boundary Conditions 

implicit real*8(a-h,o-z) 

dimension q(3,999), qn(3,999), qx(3,999), qt(3,999), 

* s(3,999), vxl(3), vxr(3), xx(999) 

c 

c nx must be an odd integer, 

nx = 101 
it = 100 
dt = 0.4d-2 
dx = 0.1d-l 
ga = 1.4d0 
rhol = l.dO 
ul — O.dO 
pi = l.dO 
rhor = 0.125d0 
ur = O.dO 
pr = O.ldO 
ia = 1 
c 

nxl = nx + 1 
nx2 = nxl/2 
hdt - dt/2.d0 
tt = hdt*dfloat(it) 
qdt = dt/4.d0 
hdx = dx/2.d0 
qdx = dx/4.d0 
dtx = dt/dx 
al = ga - l.dO 
a2 = 3. dO - ga 
a3 = a2/2.d0 
a4 = 1.5d0*al 
u21 = rhol*ul 

u31 = pl/al + 0.5d0*rhol*ul**2 
u2r = rhor*ur 

u3r = pr/al + 0.5d0*rhor*ur**2 

do 5 j = l,nx2 

q(l,j) = rhol 

q(2,j) = u21 

q(3,j) = u31 

q(l,nx2+j) = rhor 
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q(2,nx2+j) = u2r 
q(3,nx2+j) = u3r 
do 5 i = 1,3 
qx(i,j) = O.dO 
qx(i,nx2+j) = O.dO 
5 continue 

c 

open (unit=8,file= , for008 , ) 
write (8,10) tt.,it,ia,nx 
write (8,20) dt,dx,ga 
write (8,30) rhol,ul,pl 
write (8,40) rhor,ur,pr 
c 

do 400 i = l,it 
m = nx -f i - (i/2)*2 
do 100 j = l,m 
w2 = q ( 2 , j ) / q(l,j) 
w3 = q(3,j)/q(l,j) 
f21 = -a3*w2**2 
f22 = a2*w2 

f31 = al*w2**3 - ga*w2*w3 
f32 = ga*w3 - a4*w2**2 
f33 = ga*w2 
qt(ij) = -qx(2,j) 

qt ( 2 , j ) = -(f21*qx(l,j) + f22*qx(2,j) + al*qx(3,j)) 
qt ( 3 ,j ) = -(f31*qx(l,j) + f32*qx(2,j) + f33*qx(3,j)) 
s(l j) = qdx*qx(l,j) + dtx*(q(2,j) + qdt*qt(2 j)) 
s ( 2 , j ) = qdx*qx(2,j) + dtx*(f21*(q(l,j) + qdt*qt(l,j)) + 

* f22*(q(2,j) + qdt*qt(2,j)) + al*(q(3,j) + qdt*qt(3j))) 
s ( 3 , j ) = qdx*qx(3,j) + dtx*(f31*(q(l,j) + qdt*qt(l j)) + 

* f32*(q(2j) + qdt*qt(2,j)) + f33*(q(3,j) + qdt*qt(3,j))) 

100 continue 

if (i.ne.(i/2)*2) goto 150 
do 120 k = 1,3 
qx(k,nxl) = qx(k,nx) 
qn(k,l) = q(k,l) 
qn(k,nxl) = q(k,nx) 

120 continue 

150 jl = 1 - i + (i/2)*2 

mm = m - 1 
do 200 j = l,mm 
do 200 k = 1,3 
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qn(k,j+jl) = 0.5d0*(q(k,j) + q(k,j+l) + s(k,j) - s(k,j+l)) 
vxl(k) = (qn(k,j+jl) - q(k,j) - hdt*qt(kj))/hdx 
vxr(k) = (q(k,j+l) + hdt*qt(k,j+l) - qn(k,j+jl))/hdx 
qx(k,j+jl) = (vxl(k)*(dabs(vxr(k)))**ia + vxr(k)*(dabs(vxl(k))) 

* **ia)/((dabs(vxl(k)))**ia + (dabs(vxr(k)))**ia + l.d-60) 

200 continue 

m = nxl - i + (i/2)*2 
do 300 j = l,m 
do 300 k = 1,3 
q(k,j) = qn(k,j) 

300 continue 

400 continue 

c 

m = nxl -it + (it/2)*2 
mm — m - 1 

xx(l) = -0.5d0*dx*dfloat(mm) 
do 500 j = l,mm 
xx(j+l) = xx(j) + dx 
500 continue 

do 600 j = l,m 
x = q(2,j)/q(l,j) 

y = al*(q(3,j) - 0.5d0*x**2*q(l,j)) 
z = x/dsqrt(ga*y/q(l j)) 
write (8,50) xx(j),q(l,j),x,y,z 
600 continue 

c 

close (unit =8) 

10 format^ t = ',gl4.7,' it = ',i4,' ia = ',i4,' nx = ',i4) 

20 format(' dt = ',gl4.7,' dx = ',gl4.7,' gamma = ',gl4.7) 

30 format(' rhol = ',gl4.7,' ul = ',gl4.7,' pi = ',gl4.7) 

40 format(' rhor = ',gl4.7,' ur = ',gl4.7,' pr = ',gl4.7) 

50 format(' x =',f8.4,' rho =',f8.4,' u =',f8.4,' p =',f8.4, 

* ' M =',f8.4) 

stop 

end 
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Appendix B. Proof for Eq. (4.51) 


To proceed, first we shall evaluate the flux leaving each of the six quadrilaterals that 
form the boundary of a CE (see Figs. 10(a) and 11(a)). As a preliminary, note that, in 
Fig. 10(a), ^ 

area of ABGF = area of CDGB = area of EFGD = — — (5.1) 

o 

In Fig. 11(a), we have 

area of BCG A = area of DEGC = area of F AGE = ^ (5.2) 

Equations (B.l) and (B.2) can be proved easily using the information provided in Fig. 12(a). 
Moreover, because u*(x, y, t; j, k, n) is linear in x, y, and t (see Eq. (4.10)), its average 
value over any quadrilateral is equal to its value at the geometric center (centroid) of the 
quadrilateral. With the above preparations, flux evaluation can be carried out easily using 
Eqs. (4.6a)-(4.6c), (4.8), (4.10), (B.l), and (B.2). 

For each quadrilateral, the result of flux evaluation is a formula involving a x , a y , uj k , 
(u x )™ k , and {u y )J k . It can be converted to another formula involving v c , v v , uj k , («^")" fc , 
and (■«+)"*.. To carry out the above conversion, note that Eqs. (4.19), (4.20), (4.22), (4.23), 
(4.27), and (4.28) imply that 




and, for any ( j,k,n ) 6 


(5.3) 



3 

w 


( 


V 





v v 




(BA) 


Let (u x )J k , (u y )J k , . . be abbreviated as u x , u y ,..., respectively. Then Eqs. (B.3) and 
(B.4) imply that 

<■„ = ^ ) (B- 5 ) 


ha x + 

(?-*) 

3 & 
II 

+ % v v) 

(B.6) 

h(L x 

(!+») 

4 wh 
ay ~ 9a t 

(2^c T u v) 

(B.l) 
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At 

— I 

4 


u x 

-K-t 

+ “j) 


(B. 8) 

3 

wh 

{w - 6)u+ 

- (w + 

H“\-p 

(B.9) 

&X ^ X 

-f dyUy ) = 

I/^U^" + 


(«■ 10) 

W ^ 

+ 

i h 

1 Ux 2 Uj 

- = 2< 


(B-ll) 

~ 6 j 

1 h 
) u x + 2 

=r - 

r - U V 

- 2u+ 

(fl.12) 


and 


The conversion referred to above can be carried out using Eqs. (B.5)-(B.12). 

Consider Fig. 10(a). The results of flux evaluation involving the quadrilaterals that 
form the boundaries of CE r (j, fc,n), r — 1,2,3, and (j, fc,n) G fti are given here: 

(1) The flux leaving CEi (j,fc,n) through G* F f A* B* is 


2wh 




3 v s ■' / ),fc 

(2) The flux leaving CEi (j,k,n) through G'GFF' is 

2 wh / 


+ 2 I/,,) u + 2u^ - + (yc,u^ + v n u + ) 


j,k 


(3) The flux leaving CEi(_ 7 , k,n) through G'B'BG is 

2wh 


(2u( + Vr^) u-u+ + 2 + (yc u c + u *i u v') 


- 1 j,k 


(4) The flux leaving CEi(j, fc,n) through AFGB is 

(«-«C ■ 


2 wh / , , \ n “ 1 / 2 

3 V s / j + l/3,fe+l/3 


(5) The flux leaving CEi(j, k,n) through ABB 1 A 1 is 

2wh 


(t'C + 21 /,,) u - 2 + u+ - (y(W£ + !/,«+) 


n- 1/2 

j+l/3,fe+l/3 
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( 6 ) The flux leaving CEi(j, k,n) through AA' F' F is 


2u ( + u v ) U + U+ - 2 u+ - (u^U^ + Ur, 


n- 1/2 


j+l/3,fe+l/3 


(7) The flux leaving CE 2 (j,k,n) through G'B'C'D' is 


( u - 2 “< +<) iil 


( 8 ) The flux leaving CE 2 (j,k,n) through G'GBB 1 is 


( ' 2v( i 4 - u v ^ u — u~£ + 2u+ + + 


v v u v 


(9) The flux leaving CE 2 (j,fc,n) through G' D 1 DG is 


(*>< - v v ) u ~ u t ~ u v + (*W + v ’i u v ) 


(10) The flux leaving GE 2 (j,k,n) through CBGD is 


(u + 2u^ — u 1,j. 


n- 1/2 


-2/3,Jfc+l/3 


( 11 ) The flux leaving CE 2 (j,k,n) through CDD'C' is 


( 2 U( + Uqj u + u+ - 2 - (u^ 


+ + U v u t 


n-1/2 


(12) The flux leaving CE 2 (j,A:,n) through CC'B'B is 


(u n - u + u+ + u+ - (u c u+ + v v u+ 


n-1/2 


j-2/3,k+l/3 


(13) The flux leaving CE 3 (j,fc,n) through G' D' E' F' is 


(ti + «+-2«+)^ 


(14) The flux leaving CE 3 (j,fe,n) through G'GDD' is 


u — ut — ut 


J - u v + ( "C"? + u v u t 
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(15) The flux leaving CE 3 (j, k,n) through G'F'FG is 


2wh 

~9~ 


(i'C + 2 Vr^j u + 2u~£ - 




(16) The flux leaving CE 3 (jr, k,n) through EDGF is 

2wh 


u — -f 


) T2 

, 


n - 1/2 

J + l /3,fc— 2/3 


(17) The flux leaving CE 3 (j, k,n) through EFF'E' is 
2wh 


\ u < ~ v v) u + u t + u n ~ ( u < u ( + *V*j) 


1 n - 1/2 
j+1/3,*— 2/3 


(18) The flux leaving CE 3 (ji,fc,n) through EE'D'D is 

2wh 


i u < + 2l> v) u - + "v u v) 


1 n — 1 /2 
j+ 1/3, fc — 2/3 


Consider Fig. 11(a). The results of flux evaluation involving the quadrilaterals that 
form the boundaries of CE r (j,fc,n), r = 1,2,3, and (j,k,n) £ fl 2 , are given here: 

(19) The flux leaving CEi(./,fc,n) through G'C'D'E' is 


2wh 




(20) The flux leaving CEi(j,A;,n) through G'GCC 1 is 

2 wh 


u < + 2v v) u - 2u t + u v + + u v u v) 




(21) The flux leaving CEi (j,Ar,n) through G'E'EG is 

2 wh 


/ \ “I 71 

( 2l/ C + u + - 2u+ + 

*■ -I 


(22) The flux leaving CEi(j,fc,n) through DCGE is 

2 wh / , _L\ n-1 / 2 


m + «7 + 

V 4 V /j-l/3,Jk-l/3 
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(23) The flux leaving CEi(j,A:,n) through DEE'D' is 

+ 2Uv ) u + 2u t ~ u v ~ (?< u t + u i v v) 


(24) The flux leaving CEi(j, k,n) through DD'C'C is 

(2V{ + Vr^j u - + 2 - ( vc u < + u v u t) 


2wh 

~9 


n — 1/2 

j — 1 /3,fc — 1 /3 
n — 1/2 

/ — l/3,fc — 1/3 


(25) The flux leaving CE 2 (j,k,n) through G'E'F'A' is 

2 wh ( . + \ n 

— (“ +2 “< ~ 

(26) The flux leaving CE 2 (j, k,n) through G'GEE* is 

(2^ + 1/,) U + U+ - 2u+ + + u n u + ) 

(27) The flux leaving CE 2 (j,fc,n) through G* A' AG is 


2wh 

~9 




2 wh 


r -I n 

(u v -U^j U + U+ + U+ + + lSr,U+^J 


(28) The flux leaving CE 2 (j,k,n) through FEGA is 

2 wh / n . +\ n ~ 1/2 

( u — 2ut + -u+ ) 

3 V C n / j+ 2 / 3 ,k-l /3 

(29) The flux leaving CE 2 (j,fc,n) through FAA'F' is 

+ v .fl«-«++2uJ-(i<cuJ+i/,uj) 


( 2 ^ 


(30) The flux leaving CE 2 ( j , A: , n. ) through FF'E'E is 


2wh 


(j'C - ^ 


U - — u 


J - (^“c + *v4) 


n — 1 /2 
J j+2/3,fc-l/3 

n — 1 /2 

J j+2/3,k-l/3 


(31) The flux leaving CE 3 (j>, k,n) through G'A'B'C' is 

2wh 


u-u+ + 2 u+) 
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(32) The flux leaving CE 3 (ji, k,n) through G'GAA' is 


2wh 


("< - "n) u + + u v + (^f u c + ) 


-I j,k 


(33) The flux leaving CE 3 (j,k,n) through G'C'CG is 

2wh 


r -| n 

[ye + 2^7 7 ) u - 2 + u+ + ( + J4,U+ ) 

L 


(34) The flux leaving CE 3 (j,k,n) through BAGC is 

2 wh / , , \ n-1/2 


(«+«+— 2 »+)”_ i/s> 


fc + 2/3 


(35) The flux leaving CE 3 (j,^n) through BCC'B' is 
2wh 


u v ~ u <) u-u+ -u+ - (y(U+ + 1 s v u+^ 


1 »-l /2 

j-1/3, fc+2/3 


(36) The flux leaving CE 3 (j, A;,n) through BB'A'A is 


2wh 


(y< + u + 2 - u+ - + *V^) 


1 71 — 1/2 


/ — 1/3, fc+2/3 


With the aid of Eqs. (4.29)— (4.46) and (4.49a)— (4.50c), Eq. (4.51) is the result of 
(l)-(36) and Eq. (4.11). QED. 
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Appendix C. Proof for Eq. (6.51) 

As a preliminary, note that Eqs. (6.18), (6.21), and (6.27) can be used to obtain 


4 

frnt = — X fm,l {ft,q U q * + fi,q U 9 V ) 


*,9=1 


and 


4 

frnt = ~ X fm,t {ft, q U 1 x + fl,q U 9y) 


*,9=1 


In this appendix, we adopt the same convention stated following Eq. (6.32). I 
from Eqs. (6.29)-(6.32) that 


f X m,t\ 2 ( W ~ b W + b \ ( fm,t 

flj ^ V -h h 

( 1 


and 


3 a ^ \ l l / \ xy+ 

. / 


m = 1, 2, 3, 4 






3^ 

it? 




1 


w + b w — 6 

V v h 


u 


u 


m( 

+ 

TTtTf , 


m = 1, 2, 3, 4 


An immediate result of Eqs. (C.3) and (C.4) is 


X {fm,t u tx + fm,t U ty) = ^ X (/£t* U *C + -C* U t v ) » TO - 1,2, 3, 4 
£=1 *=1 

By using Eqs. (6.14), (6.16)-(6.21), and (C.1)-(C.5), it can be shown that 

3 

u mx = ( U m ( “t” U mr\ ) 


b W 


2 6 


h 


X 4 * I Umx ~ u m y ^ U mv U 


my — ^^mr) u 'm( 


( b w\ h i ~ + 

( 2 “ "g ) U mx + 2 — U mr) — 


^ + (f->)ft = 4 ^E(&+ 

^=1 


vz) 


Ut 


(C.1) 

(C.2) 

follows 

(C.3) 

(C.4) 

(C.5) 

(C.6) 

(C.7) 

(C.8) 

(C.9) 
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-Srt(»/£i+/Si)K-S) 




^=1 


i, - 2 “«' 


(C.12) 


Aft,+ 


16w/i 


mt V3 9(a/) 2 


E (/& + 2 C) (/£W< + fX<) (c.i3) 


+ (| + a) ft, = e { 2 & + /&) (/&<< + (c.i4) 


9(a<) 


fy ±—fv ±—fy 

J m J m2! - 1 - ^ J m< 


3 a / 5Z (•© /£!*) ± u e< ±u tv^^2 {ft* u q< + fi 


’ v+ u + 

i iq U qT) 


(C. 15) 


fy ± — fy zr — fy 

J m o J rax i A J mt 


I] (/mT< - /mt<) ± ± «*, ± S (/M “J + ft* U tv) 

1 = 1 L q=l J 


(C.16) 


Note that each of Eqs. (C.15) and (C.16) represents two equations. One corresponds to 
the upper signs; while the other, to the lower signs. 
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Next we shall evaluate the flux of h* m leaving each of the six quadrilaterals that form 
the boundary of a CE (see Figs. 10(a) and 11(a)). The evaluation procedure is similar to 
that described in Appendix B. For the current case, the key equations used are Eqs. (4.6a)- 
(4.6c), (6.15), (6.23)-(6.25), and (C.6)-(C.16). Futhermore, as will be shown shortly, the 
structures of the results obtained here a re very similar to those given in Appendix B. 

Consider Fig. 10(a). The results of flux evaluation involving the quadrilaterals that 
form the boundaries of CE P (j,fc,»), r = 1,2,3, and (j,k,n) € fti, are given here: 

(1) The flux of h * m leaving CEi (j,k,n) through G' F' A' B' is 


2 wh 

~1T 



+ U 


+ 

m( 


+ u 


+ 


mri J 


n 


(2) The flux of h * m leaving CEj (j,fc,n) through G'GFF' is 


2 wh 


E (At< + 2 C<) »« + 2u i - < + E (/<>* + /’>«) 

1 = 1 L 9=1 


jA 


(3) The flux of h * m leaving CEi(j,k,n) through G' B' BG is 


2wh 


E + C) 


i= 1 


U l + 


2 <4, + E (/fW< + /?>«) 


9=1 


j,* 


(4) The flux of leaving CEi (j,k,n) through AFGB is 

2 wh ( ,i \ n_1 / 2 

3~ \ Um " i+l/3,fc+l/3 

(5) The flux of h leaving CEi(j ; ,fc,n) through ABB' A' is 

^ I E (ta + 2 /s) [”' - - E (/?>£ + #«i) 

w=i L 9=1 

(6) The flux of leaving CEi(j,fc,n) through AA' F' F is 

^ ( E ( 2 C< + /a) [«« + «£ - 2 < - E («<< + tf-i) 

l*=l L 9=1 

(7) The flux of h* m leaving CE 2 (j, &,n) through G'B'C'D' is 


n- 1/2 

j+l/3,*+l/3 
n- 1/2 

► 

i+l/3,fc+l/3 


2wh 


2u^+u rn7J ^^ 
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( 8 ) The flux of h ^ leaving CE 2 (j,k,n) through G'GBB' is 
4 f 4 


2wh 


.e=i 


X ( 2 /it/ + /it) Ut U ^C + 2 “<i + X/ (//.t^tc + flq u tv; 


9=1 


(9) The flux of leaving CE 2 (j, fc,n) through G' D' DG is 


2wh 


E (/Si - C) 


k ^=1 


- «+ - < + X! <„) 


9= 1 


(10) The flux of h ^ leaving CE 2 (j,fc,n) through (75GZ? is 


2wh 


/ \ n — l/2 

I"™ + 2 "™< - "”-h_ J/s . 


3 v j — 2 /3,fc+l /3 

( 11 ) The flux of h ^ leaving CE 2 (j, fc,n) through CDD'C 1 is 


2wh 


{ X ( 2 /it/ + fi) U( + 2 “^ X (/l* U q< + 

w=i L 9=1 


( 12 ) The flux of h ^ leaving CE 2 (j,&,n) through GG'5'i? is 


2wh 


,1=1 


X (/it* - /it<) + u i< + u t - X (/£“£ + u tv) 


9=1 


(13) The flux of h* m leaving CE 3 (j,k,n) through G'D'E'F' is 

2wh 


(u m -f u 2 tt„J 

v 7 Ji* 


(14) The flux of h ^ leaving CE 3 (ji,fc,n) through G'GDD' is 


2wh 


E (/S - /&) 


^=1 


- uj - n+ + J] + /tt u tO 

9=1 


(15) The flux of h ^ leaving CE 3 (j,fc,n) through G'F'FG is 


( X (/S + 2 /ii) u * + 2 < - u t + X (//>Jc + /# 

W=1 L q=l 



- 1/2 

2/3,fc-f 1/3 
- 1/2 

-2/3,fc+l/3 
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(16) The flux of h* m leaving CE 3 (j,k,n) through EDGF is 

2 wh ( , . . \ n_1 

X 171 ~ U mC, + 2u mT) J ' 


/+l/3,fc — 2/3 


(17) The flux of h* m leaving CE 3 (j,k,n) through EFF'E 1 is 

¥ ( t (/& - e) [•* + «« + < - 1 (#•£ + fS<)] } 

lf=l L 9=1 -I > 


n- 1/2 


(18) The flux of h* m leaving CE 3 (j,k,n) through EE’D'D is 


2wh 


J ' j+l/3,*-2/3 


« r 4 

E(ta+*fls) «<-2<+«t^E(/s<c+/?> a + 0 

£=1 L q=l -I ^ j+l/3,fc-2/3 


Consider Fig. 11(a). The results of flux evaluation involving the quadrilaterals that 
form the boundaries of CE P (j,fc,n), r = 1,2,3, and (j,fc,n) £ ft 2 , are given here: 

(19) The flux of h ^ leaving CEi {j,k,n) through G'C'D'E 1 is 


2wh 


(' Um u ™< U7nr, )j,k 


(20) The flux of h * m leaving CEi(j,fc,n) through G'GCC' is 


2wh 


I X] + 2 ^m"!<) U( H + U t v + ^ (/f,t U q< + f Zq U tv) 

w=i L 9=1 


J ' j.fc 


(21) The flux of h ^ leaving CEi (j,k,n) through G' E' EG is 


2 wh 


( t, ( 2 ^< + &) [“' + < - 2 < + E (/fWc + !?,< 

I f=l L 9=1 


J ' j,k 


(22) The flux of h leaving CEi(j,fc,ri) through DCGE is 

2 wh / . . \ n— i/ 2 


U m+ U tc+ U t V )._ 1/3ik _ j/3 


(23) The flux of h* m leaving CEj (j,k,n) through DEE'D' is 
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n- 1/2 


j-l/3,fc-l/3 


n — 1/2 


j-l/3,fc-l/3 


{ 4 r 4 1 1 n 1 ^ 2 

E (/£< + 2 C<) “< + 2 < - »/, - E (/£<< + /?>«) 

£=1 L S'- 1 -I ' j—l/3,k 

(24) The flux of h ^ leaving CEi(j,k,n) through DD'C'C is 

~^\j 2 { 2 & + f^e) u t -u++ +/;>+) ] 

^ t=1 L 9=1 -J ' j-l/3,fc 

(25) The flux of h, ^ leaving CE 2 (j,k,n) through G'E'F'A' is 

2 wh ( . , \ n 

~3~ \ Um + 2Um <~ Um Vj, k 

(26) The flux of h ^ leaving CE 2 (j, fc,n) through G'GEE' is 

9~ { X + -G) U * + U l< ~ 2u X + X + fVq U tn) } 

l 1=1 L 9=1 J ) 


(27) The flux of h* m leaving CE 2 (j,k,n) through G'A'AG is 


Ut + u tc + + X {^ u t< + feJi u tn) } 

l^=i L 9=1 J ) 


(28) The flux of leaving CE 2 (j», k,n) through FEGA is 


(u m 2 


■m£ ' TTlTj 


n-1/2 

J+2/3,fc — 1/3 


(29) The flux of h* m leaving CE 2 (j,k,n) through FAA'F' is 

^ (X ( 2 /& + f£i) u * - u i + 2 u t v - X { f S u tc + ft, 

v 1=1 L 9=1 


n+ u + 
l,q u qv 


(30) The flux of h ^ leaving CE 2 (j,k,n) through FF'E'E is 

{ X (f%t - f^e) ut - u i - u t~Y, (fe« u t< + OJ?) } 

v 1=1 L 9=1 J ) 


j+2/3,* — 1/3 


n-1/2 


j+2/3,fc-l/3 
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(31) The flux of h* m leaving CE 3 (j,k,n) through G'A'B'C' is 


(^ u m 


(32) The flux of h* m leaving CE 3 (j, k,n) through G'GAA 1 is 


4 i \ n- 

(£(#-«) -« + *4 + -£ + E + 

w=i L « =l -I ' j,* 


(33) The flux of h* m leaving CE 3 (j,fc,n) through G'C'CG is 


( E (C + “< - w c + «?, + E (/& < + O 

lt=l L 9 =1 


(34) The flux of h* m leaving CE 3 (j,fc,n) through BAGC is 


( u m "+■ u m( < ^ u mi 


n — 1/2 

j-l/3,fc+2/3 


(35) The flux of h* m leaving CE 3 (j,A:,n) through BCC'B' is 

{ E (/£ - O [«' - - -4 - E (/£«£ + c 

w=i L 9=1 

(36) The flux of leaving CE 3 (j, fc,n) through BB'A'A is 

t* ( E (C + V£i) + H -<,-E (/£»£ + /’• 

I i=\ L 9 =1 


n 1 j 2 


j-l/3,fc+2/3 


U 9V 


n — 1 / 2 


j-l/3,fc+2/3 


With the aid of Eqs. (6.33)-(6.50) and (4.49a)-(4.50c), Eq (6.51) is the result of 
(l)-(36) and Eq. (6.28). QED. 
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Appendix D. Supplementary Notes 


D.l. A Discussion of Eq. (4.75) 

Here we shall prove an assertion made in Sec. 7 about the 2D a scheme, i.e., the 
backward characteristic projection of a mesh point 6 0 at the (n — l)th time level 

is in the interior of the numerical domain of dependence of the same mesh point if and 
only if Eq. (4.75) is satisfied (see Fig. 22). For simplicity, hereafter the above mesh point 
will be referred to as point 0 (not shown). In Fig. 22, the spatial projection of point 0 
at the (n — 1 )th time level is represented by point O'; while the backward characteristic 
projection of point 0 at the ( n — l)th time level is represented by point P. Without any 
loss of generality, we shall assume that j = k = 0. Thus (i) 

C = V = 0, and t — nAt (D.l) 

for point O, and (ii) 

£ — 77 — 0, and t = (n — 1 )a£ (Z>. 2) 


for point O'. 

To simplify the discussion, Eq. (4.1) is converted to an equivalent form in which (, 77, 
and t are the independent variables, i.e., 


du du du 

ai + B «ac +a ’^ 


= 0 


( 0 . 3 ) 


Here and a v are defined in Eq. (4.22). The characteristics of Eq. (D.3) are the family 
of straight lines defined by 


C = a^t + cj, and 77 = a v t + c 2 (Z>.4) 

where ci and c 2 are constant along a characteristic, and vary from one characteristic to 
another. Because points O and P share the same characteristic line, Eqs. (D.l) and (D.4) 
imply that 

C ~ — 7] = —a v At, and t = (n — 1)a< (D.5) 

for point P. Note that the temporal coordinate, i.e., t — {n — 1)a/, of points O' and P are 
suppressed in Fig. 22. 

According to the definition given in Sec. 7, the numerical domain of dependence of 
point O at the (n — l)th time level is the hexagon depicted in Fig. 22. Here the term 
‘hexagon’ refers to both the boundary and the interior. The coordinates (C? 7 /) of the 
vertices A, P, ( 7 , D, P, and F are given in the same figure. The six edges of the hexagon 
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and their equations on the (- 77 plane are 


AB 
UE 
BC 
EF 
CD 
FA 

Here the normalized coordinates £+ and 77 + are defined by 

£+ £/a£, and ri +d =rj/Ar) ( D.l ) 

As a result of Eq. (D. 6 ), a point ((,7]) is in the interior of the hexagon ABCDEF if and 
only if 

|<+ + 77 +I < 1, |77 + 1 < 1, and |<+|<1 (-D- 8 ) 

Equations (D.5), (D.7) and (D. 8 ) coupled with Eqs. (4.27) imply point P is in the interior 
of the hexagon ABCDEF if and only if Eq. (4.75) is satisfied. QED. 

D.2. The Local Euler CFL Number 

The definition of the local Euler CFL number at the point 0 (the same point defined 
in Sec. D.l) is given here. 

To proceed, consider Fig. 23. In this figure, point O' and the hexagon ABCDEF are 
also those defined in Sec. D.l. Let it, v and c be the ^-velocity, the y-velocity and the 
sonic speed at point 0, respectively. Let e x and e y be the unit vectors in the x- and the 
y- directions, respectively. Let q denote the velocity vector at point O, i.e., 

q d = ue x + ve y (-0-9) 

Let the point P depicted in Fig. 23 be at the (n - l)th time level with its spatial position 
defined by 

CLP = -qAt (D.10) 

Point P is the center of the circle depicted in Fig. 23. This circle lies at the (n - l)th time 
level and has a radius of cAt. Furthermore, it is the intersection of (i) the Mach cone [62, 

p.425] with point 0 being its vertex, and (ii) the plane with t = (to - 1)a t. For the Euler 

equations Eq. (6.10), and in the limit of At — *■ 0, this circle is the domain of dependence 
of point O at the ( n - 1 )th time level. Here a circle refers to both its circumference and 
interior. The local Euler CFL number v e at point O will be defined such that u e < 1 if 
and only if the domain of dependence of the Euler equations (i.e., the circle) lies in the 


( + + r] + = 1 
C + +V + = -i 
77 + = 1 
77 + = -l 

C + = -! 

C+ = 1 


(D.6) 


N AS A/TM— 1 998-208843 


87 



interior of the numerical domain of dependence (i.e., the hexagon ABC DEF). In other 
words, v e < 1 if and only if the normalized coordinates (C + >* 7 + ) of every point on the 
circumference of the circle satisfy Eq. (D. 8 ). 

As a preliminary, let (i) dC denote the set of the points on the circumference of the 
circle defined above, and (ii) S e denote the set of the unit vectors on the x-y plane. Then, 
for any point R £ dC (see Fig. 23), there exists an e £ 5 e such that 



PR = cAt e 

(All) 

Combining Eqs. 

(D.10) and (D.ll), one has 



O' R = (ce — q)Ai 

(D. 12 ) 

To proceed further, note that Eqs. (4.18), (4.20) and (D.7) imply that 



__ >4 _ 1 , _ w b ^ , 

vc ~ 2w (e * h Cy) 

(D.13) 

and 

_ 1 1 /^ w — b ^ 

Vrl ~2w (e * + k e < } 

(D.14) 

Let (i) £ + (C>'), £ + (P) and ( + (R) denote the values of £+ at points O', P and R, respec- 
tively, and (ii) t] + (0'), r) + (P) and rj + (R) denote the values of 77 + at points O', P and 
R, respectively. Then, because ( + (0') = r] + (0') = 0 and the gradient vectors given in 
Eqs. (D.13) and (D.14) are constant, Eqs. (D.10) and (D.12)-(D.14) imply that 


( + (P) - -At q • V<+ = (u - W + b v) 

2w h 

(D.15) 


V + (P)--Atq-V V + = -~(u+ W ~ i v) 

2 w h 

(D.16) 


C + (R) = At (ce- q) • V( + = < + (P) + cAt e • V<+ 

(D.17) 


t) + (R) = A t (ce — q) • Vr] + = T 7 + (P) -f cAt e ■ V 7 / + 

(D.18) 

and 

C + (R) + V + (R) = C + (P) + V + (P ) + CAt e • V«+ + 7 7 + ) 

(D.19) 

Note that point R is a function of e € S e . In the following, we shall evaluate the 
maxima and minima of C+(P), t] + (R) and ((+ (R) + rj+ (R)) over the range S e . To proceed, 
let 

*4 1} d = (-9 ■ ■ VC + ± c|VC + |) At (£>.20) 


1 4 2) = f (-9 • V 77 + ± c|V 77 + |) At 

(D. 21 ) 
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(£>. 22 ) 


*4 3> =' [-9 ■ V(C + + V + ) ± c|V(C + + r> + )} A t 


and 


_ def V( + 

ei = 


IVC+I’ 


- def VfJ + 

' 2 _ |V,+ |' 


and 


^3 


def V(( + + *7 + ) 

|v(C+ + ^+)| 


(£>.23) 


With the aid of Eqs. (D.13)-(D.16), (4.14) and (4.15), Eqs. (D.20)-(D.22) imply that 



■4’ - 2wh \ hu (“ + b > T “»jJ - < + (P) ± 2mh 

(£>.24) 


■4 21 - 2 Wh [hu + (W *)» T “<1 - 1 + ( p ) ± 2 wh 

(£>.25) 

and 

<4* ) - 2 wh l2Hu 2 ^=Fc i r]-C + (m’l + (J > )± 2wh 

(£>.26) 

where a£, 

at] and 



at 2 \/b 2 + h 2 

(£>.27) 


respectively, are the lengths of the three sides DF, BD , and FB of the triangle A BDF 
depicted in Figs. 12(a)-(c). Furthermore, as a result of Eq. (D.23), (i) ej is normal to 
any straight line along which is a constant, (ii) e* 2 is normal to any straight line along 
which r/~^ is a constant, and (iii) e 3 is normal to any straight line along which -f- 77 
is a constant. It follows from the above observations and Eq. (D.6) that ej, e 2 and e 3 , 

respectively, point in the directions of O' I , O' J and O' K (see Fig. 23). 

With the aid of Eqs. (D.20)-(D.23), it is easy to conclude from Eqs. (D.17)-(D.19) 

that: 

(a) For all e£ S e5 

i/W > ( + (R) > u™ (D. 28) 

with the understanding that the upper bound 1 /^. ^ and the lower bound u_ , respec- 
tively, are attained when e — t\ and e = — ej. 

(b) For all e G S e , 

v ( + ] > V + (R) > ^ ( D .29) 

with the understanding that the upper bound and the lower bound respec- 

tively, are attained when e = e 2 and e = — e* 2 . 

(c) For al 1 e 6 5 e , 

* 4 3) > C + ( R ) + v + ( R ) > 

with the understanding that the upper bound and the lower bound respec- 

tively, are attained when e = e 3 and e = — e 3 . 
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Let 


max{|j/^|, t = 1,2,3 (0. 31) 

Then Eqs. (D.24)-(D.26) imply that 

= - — - [| hu — (w + b)v | + cat]} (0. 32) 

2wh 

^ (2) = l\hu + (^ - b)v\ + ca(] (0. 33) 

2wn 

and 

“ x~~r [2 | hu — bv\ + cat ] (0.34) 

Let i/ e , the local Euler CFL number at point 0, be defined by 

v e =ma x{i/ (1) ,i/< 2 \i/ (3) } (0.35) 

Then the conclusions given in (a)-(c) coupled with Eq. (D.8) imply that the circle depicted 
in Fig. 23 lies entirely in the interior of the hexagon ABC DEF (i.e. , the analytical domain 
of dependence of point 0 lies within its numerical domain of dependence) if and only if 

v e < 1 (0.36) 


The mesh with b — 0 is used in [3]. For this special case, we have 

a( = at] — \J w 2 -f /i 2 , and at — 2h if b — 0 
As a result, Eqs. (D.32)-(D.35) imply that 
(c 4 - |u|)At At 


(0.37) 


v e — max 


{ 


w 


’ 2 wh 


[h |u| + w |v| + \/w 2 + h 2 c]| if b = 0 


(D.38) 


Note that the second component within the parentheses in Eq. (D.38) is a simplified form 
of the expression given on the extreme right side of Eq. (D.8) in [9]. As a result, v e given 
in Eq. (D.38) is identical to that given in Eq. (D.9) in [9]. 

D.3. An Existence Theorem 

Here we shall prove the following theorem. 

Theorem. At any mesh point ( j,k,n ) € fi, existence of 


y0) + 

2j (1 


and 



l = 1,2,3 


NASA/TM — 1 998-208843 


90 



is assured if the local CFL number 


u € < 2/3 


(D. 39) 


Proof: As a preliminary, we shall discuss the eigenvalues of the matrix 

M(k x ,k y ) d = k x F x + k y F* (DAO) 

Here (i) k x and k y are arbitrary real numbers, and (ii) F x and F y are the matrices formed 
by f^ t and t (see Eq. (6.13)), m,t = 1,2, 3, 4, respectively. By using (i) Eqs. (1.1), 
(1.2), (2.1) and’ (4.1)-(4.3) in [63], and (ii) the fact that two similar matrices have the 
same eigenvalues, counting multiplicity [54, p .45] , one concludes that the eigenvalues of 
M{k x ,k y ) are Aq, Aq, A+ and A_ with 


Ao 


def 


k x u + k y v 


(0.41) 


and 

A± =' A„ ± Cyjkl + (0.42) 

Note that it is assumed here that the flow variables are evaluated at the mesh point (j,k,n) 
(i.e., the point 0 referred to earlier in this appendix). 

Because F < > + and F v+ , respectively, are the matrices formed by and m,£ = 

1,2, 3,4, Eqs. (6.29), (6.31) and (4.20) imply that 


and 


pC+ _ 

4w 


pv+ _ 

4 w 




(DAO) 

(DA4) 


F <+ + pv+ 



(£>.45) 


With the aid of Eqs. (4.14), (4.15), (D.27), (D.40)-(D.45), one arrives at the following 
conclusions: 

(a) The eigenvalues of F are Aq*\ A^, A^ and A 1?^ where 


^(l) def 3A t 

0 = 4uT 



w + 

”7T 



(D.46) 


and 


x(D def x (i) 


T 


3cAtATj 
4 wh 


(DA7) 
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(b) The eigenvalues of F v+ 

are Aq 2 \ Aq 2 \ A^ and A^ where 



»(2) def 3 a* / w — b \ 

A » = 4» (“ + k V 

(£>.48) 

and 

,(2) def ,(2) _ ScAtAC 

A± - A » * to/. 

(£>.49) 


(c) The eigenvalues of (E^ + + F v+ ) are Aq 1 ^ + Aq 2 \ A^ + Aq 2 \ A^ and A^ where 


\ ( 3 ) def .(1) \ (2) 3CA <AT 

A± — A 0 + A 0 =F 4wh 


{D. 50) 


Let (i) Aj, A 2 , X n be the eigenvalues of any nx n matrix A, and (ii) I be the 
n x n identity matrix. Then the eigenvalues of the matrix I — A are 1 — Aj, 1 — A 2 , . . 

1 — A„. As a result, Eqs. (6.33), (6.36), (6.39), (6.42), (6.45) and (6.48) coupled with the 
above results (a)-(c) imply that: 

(d) The eigenvalues of Sj^ + are 1 — A^ — A^ 2 *, 1 — A^ — A^, 1 — A^ and 1 — A^ while 
the eigenvalues of Ej 2)+ are 1 + A^ + A^ 2 *, 1 + A^ + A^ 2 , 1 + A ^ and 1 + 

(e) The eigenvalues of arel + A^, 1 + Aq 1 ^, 1 + A^ and 1 + A^\ while the eigenvalues 

of E^ 2)+ are 1 - aJ 1} , 1 - A^, 1 - A ( + J) and 1 - A ( _ 1) . 

(f) The eigenvalues of arel + AQ 2 \ 1 + Ag 2 \ l + A^ and 1 + A^, while the eigenvalues 

of Sj 2)+ are 1 - X ( Q 2 \ 1 - \ {2 \ 1 - A ( + 2) and 1 - A ( _ 2) . 

Note that the matrices referred to in (d)-(f) are nonsingular, and therefore their inverses 
exist, if the eigenvalues of these matrices are nonzero [54, p.14]. To complete the proof, 
we need only to show that these eigenvalues are nonzero if v € < 2/3. 

To proceed, note that, because c > 0, it follows from Eqs. (D.24)-(D.26) that 

> C + (P) > V- \ and > r/ + (P) > (D. 51) 

and 

v { + ] > <+(P) + V + (P) > ^- 3) (£>.52) 

With the aid of Eqs. (D.31), (D.35), (D.51) and (D.52), Eq. (D.39), which is equivalent to 
(3/2)i/ e < 1, implies that 

| i^l < 1, i — 1,2,3 (£>.53) 

and 

^|C + (P)|<1, ^|7 ? + (P)|<1, and ^|< + (P) + 7 7 +(P)|<l (£>.54) 
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Next note that Eqs. (D.15), (D.16), (D.24)-(D.26) and Eqs. (D.46)-(D.50) imply that 


A 


(0 _ 
± - 



£ = 1,2,3 


(D. 55) 


and 

A<'U-?C + (n 4 J) = -?l, + (n and Aj, 1 ' + A^ 2 > = _|(C+(J>) + »? + (-P)) (X>.56) 

It now follows from Eqs. (D.53)-(D.56) that each one of the eigenvalues listed in (d)-(f) 
has the form of 1 ± x with \x\ < 1 if u e < 2/3. Thus these eigenvalues are all positive if 
*/ e < 2/3. QED. 
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(a) Profiles at t = 0.2 


Figure 1. The t'E/SE solution of the extended Sod’s problem 
using the boundary conditions Eqs. (2.60) and (2.67) 
(A/ = 0.004. A.r = 0.01, C’FLw 0.88. < = 0.5, o = I) 
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Figure 4. — The SEs and CEs. (a) The staggered space-time mesh, (b) SE(j,n). 
(c) Alternative SE(j,n). (d) CE-(j,n). (e) CE+(j,n). (0 CE(j,n). 
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Figure 5. — A spatial domain formed from congruent triangles, showing 
the spatial projections of the mesh points. 
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A 


Figure 6. — 
positions o 


NASA/T 



Figure 7. — The relative spatial positions of the mesh points 
e ft] and the mesh points € fl2 (dash lines are 
spatial boundaries of the conservation elements 
depicted in figs 10(a) and 1 1(a)). 


(- 2 , 4 ) 



Figure 8. — The spatial mesh indices (j, k) of the mesh points e Hj 
(n = ±1/2, ±3/2, ±5/2, ••-). 
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Figure 9. — The spatial mesh indices (j, k) of the mesh points e H 2 
(n = 0, ±1, ±2, •••)■ 
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k 


CE 1 (j. k, n) = box GFABG'F'A'B' SE (j, k, n) = the union of four 

CE 2 (j, k, n) = box GBCDG'B'C'D' planes A'B'C'D'E'F', GBB "G ", 

CE 3 (j, k, n) = box GDEFG'D'E'F' GDD "G ", and GG"F"F and their 

immediate neighborhoods. 

G' = (j, k, n) e 

A' = (i + -|,k + -l,n), B=0-y,k + |,n). C' = 0 - |, k + 1, n), 

D=(i_ i’ k "3’ n) ’ E=(j+ 3 ,k "|’ n) ’ F ' = 0+|,k- J,n) 

Figure 10. — (a) Conservation elements CE r (j, k, n), r = 1, 2, 3, for any (j, k, n) 
e ft], (b) Solution element SE (j, k, n) for any (j, k, n) e ftj. 
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t 



CE-j (j, k, n) = box GCDEG'C'D'E' 
CE 2 (j, k, n) = box GEFAG'E'F'A' 
CE 3 (j, k, n) = box GABCG'A'B'C' 


t 



./ 


(b) 

SE (j, k, n) = the union of four 
planes A'B'C'D'ET', GG"A"A, 
GCC"G", and GG'E'E and their 
immediate neighborhoods. 


A' = (j + -|,k + l,n), 
D =(i-i k-i n), 


G' = 0, k, n) € n 2> 

B' = (i-i k + |,n), C' = (j-|,k + -l.n). 
E =0+ 3’ k_ 3’ n) ’ F = (j+ 3 ,k_ 3’ n) 


Figure 11. — (a) Conservation elements CEr (j> k, n), r = 1, 2, 3, for any 
(j, k, n) e fl2- (b) Solution element SE (j, k, n) for any 
(J, k, n) e ft2- 
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(0 - §)A£, (k + 1) Ann) ((j - i)At, (k + f ) Ar|) 

J o 3 3 



Figure 12. — Geometry of the hexagon ABCDEF. (a) Relative positions of the vertices 
in terms of (x, y). (b) Relative positions of the vertices in terms of (j, k). 
(c) Relative positions of the vertices in terms of (£, m). 
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Figure 13. — (a) The mesh points (j, k, n), (j + 1/3, k + 1/3, n - 1/2), 

(j - 2/3, k + 1/3, n- 1/2) and (j + 1/3, k - 2/3, n - 1/2) 
that belong to ftp (b) The mesh points (j, k, n), (j - 1/3, 
k - 1/3, n - 1/2), (j + 2/3, k - 1/3, n - 1/2), and (j - 1/3, 
k + 2/3, n - 1/2) that belong to 1^2 • 
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Figure 19. — The functions Xr(e), r = 1, 2, 3. 



A 


Figure 20. — Spatial projection of part of a 3D space-time mesh, 
showing the construction of a CE. 
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( 0 , 1 / 2 ) ( 1 / 2 , 1 / 2 ) 



(b) (0,0) (1/2,0) 

Figure 21. — Concept of dual space-time meshes, (a) The dual space-time mesh. 

(b) A rectangular space-time region shared by CE_( 1 /2, 1 /2 ) and 
CE+(0,l/2). 
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Figure 22. — The numerical and analytical domains of dependence 
associated with the 2D a-scheme. 



F (At, -At|) 


Figure 23. — The numerical and analytical domains of dependence 
associated with the 2D CE/SE Euler solvers. 
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